Papers

Let $x_1,\ldots,x_n$ be independent standard Gaussian vectors in $\mathbb{R}^d$. An ellipsoid fit is a matrix $S \succeq 0$ such that $x_i^\top S x_i =d$ for every $i$, so that all the points lie on the boundary of the centered ellipsoid $\{ x : x^\top S x = d\}$. Saunderson, Parrilo and Willsky conjectured that, as $n,d \to \infty$, this semidefinite feasibility problem undergoes a sharp transition at $n \sim d^2/4$. We prove this conjecture. If $\lim \sup n/d^2 = \alpha^* <1/4$, then, with probability tending to one, an ellipsoid fit exists; moreover, one can choose $S$ with all eigenvalues in a fixed interval $[\lambda_- , \lambda_+] \subset (0,\infty)$ depending only on $\alpha^*$. Conversely, if $\lim \inf n/d^2 > 1/4$, then, with probability tending to one, no ellipsoid fit exists, without any spectral restriction. Our proof builds on the Gaussian-equivalence framework developed by Bandeira and Maillard (2025) and closes the two gaps left open in their work: establishing exact fitting and removing the operator-norm constraint. On the satisfiable side, the new ingredients are a head-tail decomposition of the dual vector, exact correction of the sparse head constraints, and a Gaussian comparison principle for the low-influence tail. On the unsatisfiable side, we split a candidate into a low-rank spectral head and a Schatten-3 diffuse bulk, Gaussianize the bulk conditionally on the head, and apply a projected Gordon escape argument. The threshold is governed by the statistical dimension $d(d+1)/4$ of the positive semidefinite cone.

Powerful black-box predictive models have motivated many proposals for combining observed data with predictions to perform valid statistical inference. Despite this progress, the field lacks a unifying principle that explains how hypothesis tests should integrate data and predictions in a way that is both valid and efficient. In this work, we address this gap in the high-dimensional Gaussian sequence model. We characterize the information-theoretic limits of inference with black-box predictions when their accuracies are known and, for orthogonal predictions, when they are unknown. Building on these characterizations, we develop practical hypothesis tests that adapt to the unknown accuracies of the predictions while benefiting from strong alignment among them.

Many modern testing problems are sequential along two axes: new hypotheses may arrive over time, while evidence for hypotheses already under consideration continues to evolve and may be inspected at arbitrary stopping times. We develop dynamic $e$-closure for this setting. At a global stopping time the active true-null intersection is random. Future-extension coherence allows its certificate to be compared with that of a fixed terminal intersection, yielding simultaneous stopped-FDR control. If the certificates are also time-monotone, the resulting closure controls simultaneous SupFDR and is setwise persistent. Conversely, every procedure satisfying either criterion is contained in a dynamic closure generated by canonical normalized-loss processes. For pointwise mergers, fixed-dimensional admissibility is equivalent to ordinary arbitrary-dependence $e$-merging. Coherence across horizons then forces a single globally summable weight sequence and exact neutrality under padding by the $e$-value one; on a countably infinite hypothesis universe, this rules out nontrivial symmetric mergers in the admissible pointwise class. The theory extends from FDP to bounded losses that are monotone in the possible true-null configuration and local in the reported action. We also give a coherence counterexample, persistent constructions, and a globally valid shared-control model.

Uniform stability is a classical tool for controlling the generalization error of a learning algorithm. Bousquet, Klochkov, and Zhivotovskiy (2020) showed that the problem can be reduced to a moment inequality for a sum of weakly interacting functions of independent random variables. Their bound contains an additional factor $\log n$, and they asked whether this factor can be removed. We answer this upper-bound question affirmatively. More specifically, let $Z=(Z_1,\ldots,Z_n)$ have independent coordinates and let $g_i(Z)$ satisfy $\mathbb E[g_i(Z)\mid Z_{-i}]=0, \ \left| \mathbb E[g_i(Z)\mid Z_i]\right|\le M, \ \text{for every } i = 1, \dots, n, $ where $Z_{-i}$ denotes all coordinates except $Z_i$. Assume additionally that changing any coordinate $Z_j$, $j\neq i$, changes $g_i$ by at most $\beta$, we prove that, for every $p\ge2$, for every $p\ge2$, $$ \left\| \sum_{i=1}^n g_i(Z)\right\|_p \le 16pn\beta+M\sqrt{2pn}. $$ This removes the $\log n$ factor from the previous bound and matches the lower bound of Bousquet, Klochkov, and Zhivotovskiy up to universal constants in the range covered by their construction. Our proof first establishes the required estimate on the Rademacher cube, then transfers it to arbitrary product distributions by a two-copy randomization argument.

We consider the Langevin diffusion $dX_t = - \beta \nabla V(X_t) dt + \sqrt{2} dB_t$ for a general nonnegative real-analytic potential $V$ and a large parameter $\beta$. In the large-$\beta$ limit the process is confined to the zero set of $V$, assuming that it starts there. We derive a candidate limiting evolution on the zero set. To do so, the zero set is partitioned into strata according to a measure of local codimension known as the local learning coefficient and its multiplicity. It is then shown that the Dirichlet form associated with $X$ converges in a certain sense to a hierarchy of Dirichlet forms corresponding to a stochastic evolution on the strata. This evolution is strongly biased toward higher-dimensional, or "more singular", strata. This result is motivated by a question from Watanabe's singular learning theory regarding the learning dynamics of overparameterized statistical models and the generalization puzzle in deep learning. The result suggests a mechanism for the observation that stochastic gradient methods tend to be biased toward singular solutions that generalize well.

We study estimation of tail heterogeneity for non-identically distributed extreme observations subject to random right-censoring. In the uncensored setting, such heterogeneity is described by the event scedasis function, which measures the relative contribution of different design points to the upper tail. Under censoring, however, the observed tail heterogeneity is contaminated by the censoring scedasis functions, and applying uncensored techniques targets the wrong object. We propose a Beran-type estimator of the relative event scedasis, which is consistent under mild conditions. To obtain these results, survival analysis representations at an upper order statistics are extended to the non-identically distributed case; specifically, we develop conditional Nelson–Aalen and Beran theory on increasing intervals whose random endpoint is dominated, with probability tending to one, by a deterministic high local quantile. In particular, we derive a martingale array representation of the conditional Nelson–Aalen estimator with explicit error bounds depending only on the sample fraction and the bandwidth. Simulations demonstrate the finite-sample performance of the method, and an application to French property-casualty insurance claims illustrates how heterogeneous censoring can distort naive scedasis estimates.

In this paper, we study kernel ridge regression using the data-driven diffusion maps (DM) kernel, which is constructed through algebraic manipulations of the diffusion maps algorithm. Under the assumptions that the data lie on a manifold and are sampled uniformly, we prove that the appropriately scaled DM kernel converges uniformly to the heat kernel on the manifold for sufficiently large times as the dataset size increases and the kernel bandwidth is scaled appropriately. Consequently, the limiting Reproducing Kernel Hilbert Space (RKHS) induced by the DM kernel coincides with the RKHS associated with the heat kernel on the manifold. We further show that the RKHS induced by the DM kernel is isometrically isomorphic to the RKHS of the Gaussian kernel, which is continuously embedded in the RKHS of a Matérn kernel whose norm is equivalent to an appropriate Sobolev norm. This result implies that standard risk bounds for kernel ridge regression applicable to Matérn kernels also apply to the DM kernel. Finally, we provide numerical results that (1) validate the convergence of the heat kernel approximation, (2) demonstrate the greater expressiveness of the DM kernel compared to the Gaussian kernel for supervised learning over a larger class of functions on manifolds with boundary, and (3) demonstrate the advantage of the DM kernel over the Gaussian kernel in learning functions with varying frequencies and co-dimensions.

In this paper, we consider the problem of estimating conditional expectations as an ill-posed inverse problem. We propose a solution based on a generalization of Vapnik's theorem \cite[Theorem 7.2]{Vapnik1998} for solving stochastic ill-posed problems in Hilbert spaces. As an application, we derive a new upper bound for sample errors of conditional expectation estimation.

How expressive is prompting a transformer? Answering this question is important for separating the roles of prompting, architecture, and pretraining in transformer models, and for determining whether task-specific behavior must be stored in model weights or can instead be induced at inference time through the prompt. We show, in an approximation-theoretic sense, that pretraining is optional: a single-layer softmax attention network with random, untrained weights can approximate any Hölder function on a compact manifold when steered by an appropriate soft prompt. Guided by the connection between softmax attention and kernel methods, we construct explicit soft prompts (a prompt per target function, independent of the query) as solutions to linear systems matching attention logits to Gaussian kernel exponents, under which the frozen transformer emulates the classical Nadaraya-Watson kernel estimator. The construction requires only a mild rank condition on the weights, which we show holds almost surely under Gaussian initialization. The prompted network inherits the theoretical guarantees of kernel regression, leading to universal approximation theorems with minimax-optimal rates that depend on the intrinsic dimension. We further quantify the cost of prompting, exposing a tradeoff between the norm of the constructed soft prompt tokens, prompt length, and hidden dimension. Numerical experiments corroborate the constructions and predicted rates.

High-dimensional kernel density estimation suffers from the curse of dimensionality. This study proposes a hybrid density estimator defined as the product of a joint density over a pre-specified subset of dimensions and marginal densities for the remaining variables, instead of estimating the full high-dimensional density directly. Under this framework, given the observed data, we select a subset of variables for joint density construction to improve estimation accuracy, while modeling the remaining variables through their respective marginal densities. This construction involves a trade-off between the approximation error induced by simplifying the dependence structure and the variance reduction achieved by lowering the dimension of the joint density. We employ a genetic algorithm to efficiently identify such variable subsets. Theoretical and numerical results demonstrate that the proposed estimator can outperform conventional full-dimensional kernel density estimation when this trade-off is appropriately balanced.

A score-driven filter multiplies its scaled log-likelihood score by a scale parameter. We call this coefficient the gain and learn it online. Given the current state and realised scaled score, each admissible gain selects a reachable next state and predictive density. A scalar gain moves along a line; diagonal gains control coordinatewise transmission and may change direction. We evaluate gain selection using a one-step predictive Kullback-Leibler objective. In the scalar unscaled case, the negative consecutive-score product is a stochastic gradient; the positive product used in accelerated recursions is a descent direction. Positive scalar score scaling changes only the effective learning rate. Strictly increasing, continuously differentiable gain links with positive derivative induce mirror-descent geometry, while persistence adds a Bregman pull towards a reference gain. Under convexity, compactness, integrability, and schedule conditions, projected and discounted mirror updates satisfy dynamic-regret bounds relative to time-varying, current-information comparators. Simulations isolate score scaling, link geometry, persistence, and coordinatewise gains. Across twelve equity indices, the bounded discounted-logistic gain records a lower out-of-sample mean negative log score than the constant gain in eleven markets, although market-level evidence is mixed. It also avoids the extreme transients of the numerically capped exponential-link benchmark.

High-dimensional Principal Component Analysis (PCA) has become an essential tool in modern data analysis, offering dimensionality reduction and feature extraction. However, the presence of missing data introduces significant challenges, distorting the performance of PCA and complicating statistical inference. In this paper, we study the asymptotic behavior of PCA under a spiked population model with missing observations, leveraging recent advances in random matrix theory. We demonstrate that while the spiked sample eigenvalues exhibit asymptotic normality, the limiting parameters differ substantially from those in the complete data case, reflecting the non?trivial influence of the missing data mechanism. As an application of our results, we propose a test to evaluate the independent structure of a spiked population.

Cells often estimate chemical concentrations from only a few diffusing molecules, whose stochastic motion limits sensing precision. We consider a perfect monitoring instrument that records the instantaneous molecular density throughout a finite region without perturbing the concentration field or distinguishing molecular identities. The standard Berg-Purcell estimator weights all positions uniformly. Here we derive the optimal spatial weighting within a large class of unbiased estimators for both concentration and gradient sensing. Variance minimization maps to an electrostatic problem. Surprisingly, although the instrument monitors the entire volume, the optimal estimator assigns all weight to its boundary, reducing the uncertainty in both concentration and gradient sensing. We extend the results to arbitrary geometries and spatial dimensions.

An expert who supplies examples of a quantity often also signals how plausible each one is; when is that signal worth using? We study eliciting a Bayesian prior from an expert who provides example points together with their approximate likelihoods. We propose fitting the prior by least squares–minimizing the squared discrepancy between a parametric density and the elicited likelihoods–which defines an M-estimator that remains well posed even for families whose moments do not exist. We establish consistency and asymptotic normality, and prove–under explicit regularity conditions, comprising a well-separation and a uniform-concentration requirement that we verify for the families considered–a non-asymptotic O(1/sqrt(n)) Berry–Esseen bound on its sampling distribution, uniform and nonuniform, by extending a result of Pinelis for maximum-likelihood estimators to the M-estimation setting. We then relax the assumptions that most limit the method in practice. Experts need not report on the density's own scale: an unknown reporting scale can be profiled out in closed form and estimated jointly, at no asymptotic cost for location families. The theory extends to multivariate parameters, where a directional Berry–Esseen bound follows from the multivariate delta method applied to a smooth implicit proxy for the estimator. An additive error floor in the noise model removes a degeneracy in the optimal design, making optimal designs interior. Simulations for normal and beta families confirm the predicted n^-1 error rate and the accuracy of the normal approximation at moderate sample sizes. Finally, we compare the estimator with the sample-only maximum-likelihood baseline, derive an explicit threshold on the expert's reporting noise below which the elicited likelihoods provably reduce estimation error, and calibrate that threshold against eleven datasets of human frequency judgments.

Kernel spectral clustering with a single bandwidth can be inadequate for data exhibiting multiple characteristic pairwise-distance scales, a problem particularly prevalent in the high-dimensional regime. We address this issue through a multi-kernel formulation that aggregates kernels with different bandwidths. The bandwidths are selected as prescribed empirical quantiles of the pairwise squared distances, thereby capturing the relevant distance scales without requiring prior population-scale information. We develop a rigorous theoretical analysis of the resulting method under a general high-dimensional, multi-scale mixture model with heterogeneous cluster centers and covariance geometries. We construct a blockwise constant, low-rank informative approximation to the empirical multi-kernel matrix and establish row-wise $\ell_{2,\infty}$ perturbation bounds for its leading spectral components, as well as for the associated normalized Laplacian matrix. These bounds yield observation-level control of the spectral embedding, which is more informative than conventional global eigenspace perturbation estimates. Under suitable eigen-gap and cluster-separation conditions, we show that approximate $K$-means applied to the multi-kernel spectral embedding achieves exact recovery with high probability.

We study the conditionally resampled sliding-window count kernel associated with the empirical counts of length-$n$ windows from a stationary finite-state reversible Markov chain. Although the resulting count process is generally not Markov, its stationary one-step conditional law defines a genuine Markov kernel. For every fixed strictly positive reversible kernel \(P\) on a finite state space, we present a Poincaré inequality for the induced count kernel $\tP_n$ of length $n$. In other words, we derive the lower bound of the spectral gap $\Gap(\tP_n)$ of $\tP_n$ as \[ \Gap(\tP_n)\ge \frac{c(P)}{n}, \] where \(c(P)>0\) depends only on \(P\). The proof combines a martingale oscillation inequality for the stationary path law with a direct comparison of coordinate oscillations to the Dirichlet form of the count kernel. A linear statistic of the count vector gives the matching \(O(1/n)\) upper bound, so for every fixed strictly positive reversible \(P\) one has \(\Gap(\tP_n)=\Theta_P(1/n)\). The resulting count-space Poincaré inequality yields a local-to-global variance bound for finite-window count statistics and, together with a general matrix-concentration principle, operator-norm concentration for matrix-valued empirical averages.

Fractional elliptic stochastic partial differential equations (SPDEs) are widely used in statistics and machine learning for computationally efficient and flexible modeling of Gaussian random fields. The computational efficiency of the SPDE approach relies on finite element approximations combined with numerical quadrature and mass lumping, which enable sparse matrix methods during inference. Although many works have studied finite element approximations of fractional SPDEs, the effect of the mass lumping and quadrature approximations used in practice has not been fully analyzed. To fill this gap, we derive convergence rates for numerical approximations of fractional SPDEs based on finite element discretizations combined with numerical quadrature and mass lumping. Specifically, we obtain explicit convergence rates for the mean-squared error of the covariance function in a general framework that covers the main settings where mass lumping is used in the SPDE approach. We also analyze non-stationary variance-control factors of the form $L^\beta(\tau u)=\mathcal{W}$, where $\tau$ is spatially varying, and derive covariance error estimates showing how the regularity of $\tau$ affects the convergence rate. As specific examples, we provide results for random fields on bounded Euclidean domains, Riemannian manifolds, and metric graphs. Numerical experiments are presented that confirm the theoretical results.

We study identifiability of mixing measures in infinite mixture models. We show that, in many common cases, lack of identifiability can be characterized in terms of certain differential structures of the kernel family with respect to its parameters. In our main results, we prove that when the kernel is annihilated by a non-trivial differential or difference-differential operator over the parameter space, there exist infinitely many distinct mixing measures yielding the same mixture density. We give verifiable conditions for such operators to exist, covering many common cases, including the location-scale Gaussian, location-scale Student-t, Gamma, Beta, Dirichlet, negative binomial and non-central Chi-squared families. Furthermore, our conditions apply to any exponential family whose parameter dimension exceeds the dimension of its sufficient statistic and, more generally, to kernels with polynomially-growing score functions. We complement our results with a minimax lower bound on the estimation error for the mixing measure in the Wasserstein distance under non-identifiability. On the flips side, we describe three classes of kernels for which identifiability is preserved and nonparametric statistical inference remains possible.

We study variable screening for high-frequency continuous-time regression, targeting the support of the spot regression coefficient. We establish uniform entrywise bounds for local estimates of the marginal correlation vector and covariate correlation matrix under coordinatewise truncation and finite-variation jumps. Under a local averaged $\alpha$-Hölder condition, the rate is $\{(\log p)/n\}^{\alpha/(2\alpha+1)}$, up to logarithmic factors. This yields sure screening for the marginal support. We also derive a list-size minimax lower bound that matches the localization and multiplicity order when the retained list size $d$ satisfies $\log(p/d)\asymp\log p$. For regression-support recovery under local dependence, we propose partial-residual spot iterative screening with pathwise Lepski adaptation. Under sparse eigenvalues and standardized beta-min, the procedure has the sure screening property and requires only correlation columns associated with the current active set. Simulation studies and a high-frequency factor application illustrate its performance.

Approximate message passing (AMP)-type algorithms are widely used for signal recovery in high-dimensional noisy linear systems. Recently, a framework called memory AMP (MAMP) was introduced, offering a new approach to incorporating memory terms within AMP algorithms. Building on this, a low-complexity gradient descent MAMP (GD-MAMP) was proposed for right-unitarily invariant matrices. In this paper, we first address an overflow problem in GD-MAMP caused by intermediate variables exceeding the floating-point range, which typically occurs when the condition number is large. Second, we propose two low-complexity variants of GD-MAMP: one replaces full-length memory with partial memory, while the other reduces the number of matrix-vector products per iteration by $1/3$ (from three to two). Neither degrades the convergence speed notably. Third, we develop a general gradient-based formulation for designing MAMP algorithms. This formulation recovers warm-started conjugate gradient VAMP (WS-CG-VAMP) as a special case. Furthermore, we show that the computation of the orthogonalization parameters in this formulation can suffer from catastrophic cancellation, which explains the finite-precision instability of WS-CG-VAMP. Finally, we derive an equivalent reformulation, termed WS-CG-VAMP(r), which reduces the number of matrix-vector products by up to $50\%$. Measured by matrix-vector products, GD-MAMP converges faster for small condition numbers, whereas WS-CG-VAMP(r) converges faster for large ones under high-precision arithmetic but may diverge in IEEE double precision due to catastrophic cancellation.

Ranking data arise in scientific and machine learning applications, including recommendation systems, information retrieval, voting, marketing, and AI preference ranking from human feedback. Existing statistical work has primarily focused on inference tasks such as preference estimation, rank aggregation, and ranking prediction. However, generating realistic synthetic rankings from an observed population is important for privacy-preserving data sharing, benchmark construction, simulation, and uncertainty quantification. This task is challenging because rankings are high-dimensional combinatorial objects with non-Euclidean dependence structures, while ranking populations often exhibit substantial preference heterogeneity. We propose a framework for population-level generative modeling through a latent preference simplex embedding. It estimates a low-dimensional latent preference simplex through a likelihood-based ranking model, leverages flow matching to learn the population distribution of latent preferences, and generates new rankings through the fitted probabilistic ranking model. We show that ranking generation admits an oracle reduction to latent distribution learning and derive finite-sample generative guarantees that clarify how the number of items, ranking length, and latent dimension affect accuracy. Experiments on synthetic and real datasets demonstrate improved population-level fidelity and provide a statistically interpretable representation of preference heterogeneity.

We study constrained statistical learning over infinite-dimensional hypothesis classes in the fully nonconvex setting, and establish universal PACC learnability of the solutions of dual algorithms: Probably Approximately Correct on Constraints, guaranteeing optimality and constraint satisfaction at once. This strengthens near-PACC results, whose feasibility residual no amount of data can remove. Optimality is caught between generalization, governed by Rademacher complexity and favoring small classes, and strong Lagrangian duality, which rests on Lyapunov convexity for vector measures and needs decomposability, a demand pulling the other way. We reconcile the two by posing the population problem over a universal RKHS $\mathcal{H}_K$, dense in a decomposable envelope, and learning over norm balls of growing radius. This yields the Tikhonov complexity $\mathfrak{T}^{\varepsilon}_{n}$, the least RKHS norm reaching an $\varepsilon$-optimal Lagrangian level set; we prove it finite, obtain exact learnability of the optimal value, and make the sample threshold explicit and polynomial in $1/\varepsilon$ under a source condition. Feasibility is harder: absent convexity the Lagrangian may not attain its infimum, and dual information pins down only an averaged constraint-risk vector, not the risks of any returned predictor. We introduce the closure-realization gap $\varepsilon^\star_\infty$, an index of how well $\mathcal{H}_K$ retrieves feasible solutions from dualization; it is a property of the problem, not of a modeling choice. Learnability is exact when $\varepsilon^\star_\infty=0$, in particular under dual differentiability, and near-PACC with residual exactly $\varepsilon^\star_\infty$ otherwise. Finally, no distribution-free threshold exists already in the unconstrained specialization, so universality is the canonical frame for dual algorithms over large hypothesis classes.

Throughout the last decade, Gaussian universality has been widely studied for high-dimensional estimation problems. Most of the literature focuses on i.i.d. sensing matrices or accounts for special forms of dependence, such as block dependence or other specific row/column dependencies. More general simultaneous row and column mixing has not yet been fully studied. In this paper, we focus on that setting. We prove a Gaussian universality theorem for the lasso in the sparse regime, where the non- Gaussian covariates have linearly dependent rows and columns. To the best of our knowledge, our setting permits a broader simultaneous row and column dependence structure than those treated in much of the prior universality literature. Numerical illustrations for various sparse profiles support the universality claims of this paper.

We study pointwise convergence rates of purely random partition estimators in nonparametric regression, where the partition – into hyper-rectangles by purely random trees, or into Voronoi cells by prototype rules – is built independently of the responses. Our analysis rests on a single geometric criterion, shape regularity, relating the diameter of a cell to its volume, which is shown by Bettinger, Portier and Saumard (2026) to be necessary and sufficient, up to logarithmic factors, for achieving the minimax rate $n^{-1/(d+2)}$. We show that centered and uniform trees are not shape-regular – their cells' aspect ratio grows exponentially with the number of splits with probability bounded away from zero – explaining the super-logarithmic corrections in their error bounds, whereas Mondrian trees, whose splits adapt to the current cell geometry, are shape-regular in probability and attain the minimax rate. The same analysis applied to Voronoi partitions yields the first pointwise concentration bounds for Proto-NN, resolving an open problem of Györfi and Weiss (2021), and shows that OptiNet achieves the minimax rate with markedly better success probability – even almost surely, for a suitable choice of parameters – thanks to its $\eta$-net construction.

This chapter deals with nonparametric inference for quantiles from a Bayesian perspective, using the Dirichlet process. The posterior distribution for quantiles is characterised, enabling also explicit formulae for posterior mean and variance. Unlike the Bayes estimator for the distribution function, our Bayes estimator for the quantile function is a smooth curve. A Bernstein–von Mises type theorem is given, exhibiting the limiting posterior distribution of the quantile process. Links to kernel-smoothed quantile estimators are provided. As a side product we develop an automatic nonparametric density estimator, free of smoothing parameters, with support exactly matching that of the data range. Nonparametric Bayes estimators are also provided for other quantile-related quantities, including the Lorenz curve and the Gini index, for Doksum's shift curve and for Parzen's comparison distribution in two-sample situations, and finally for the quantile regression function in situations with covariates.

Causal models in the same Markov equivalence class are, in the absence of strong assumptions, statistically indistinguishable at any sample size. For modest sample sizes there is also the possibility that several such classes are compatible with the data, pointing to a confidence set of causal models as the appropriate presentation of evidence. Non-identifiability of Gaussian causal models from the same Markov equivalence class corresponds to a plurality of inverse-covariance representations in the unconstrained parametrisation of Cox and Wermuth (1993). This parametrisation, avoiding conic constraints that would otherwise complicate distributional approximations, facilitates construction of, and theoretical analysis for, a confidence set of causal models based on standard likelihood theory. By drawing on a geometric formulation of Evans (2020), we provide insight into which causal models are most likely to be included in the confidence set when the true parameter values of the generating model are at the borderline of detectability, delineating those models whose inclusion probabilities are stable at the nominal level, slowly decaying with sample size, and quickly decaying with sample size. We also study settings in which two or more causal models are in operation, exploring how the mixture weights, and the geometry of the models in the mixture relative to candidate models, interact with the inclusion probabilities. The purpose of the paper is to probe, from the perspective of structural properties of the true causal mechanism, the limits of what is achievable in causal inferential settings.

We determine necessary and sufficient conditions for complete monotonicity of the Srivastava–Tomovski extension, a generalized Mittag–Leffler kernel arising in fractional calculus, while separating the defining generalized-Wright series from its positive-axis realization when the series is not entire. For positive parameters, let \(\Delta=1+\alpha-\kappa\). In the entire regime \(\Delta>0\), \(E_{\alpha,\beta}^{\gamma,\kappa}(-x)\) is completely monotone on \((0,\infty)\) if and only if \(\alpha\leq\kappa\) and \(\kappa\beta\geq\alpha\gamma\). For arbitrary positive parameters, the corresponding positive-axis realization satisfies the same criterion; when \(0<\alpha<\kappa\), it is defined by a second-kind Wright kernel and remains meaningful through the finite-radius and zero-radius phases of the defining series. We also identify the Bernstein measure in every admissible case. In the interior region it is, after normalization, the pushforward under \(t\mapsto t^\kappa\) of a power-biased Wright distribution; at \(\alpha=\kappa\) it becomes a powered beta law for \(\beta>\gamma\) and the atom \(\Gamma(\gamma)^{-1}\delta_1\) for \(\beta=\gamma\). Necessity in the excluded regions follows from a Mellin-transform zero, uniqueness of finite signed Laplace transforms, and moment-support asymptotics. The Prabhakar criterion is recovered by setting \(\kappa=1\).

In the classical normal means problem, independent train–test folds can be constructed by perturbing the data with normal randomization. Averaging over $K$ such folds yields a cross-validation estimator whose bias depends on the marginal distribution of the randomization variables, while its variance depends on their joint distribution. This raises the questions: which joint law is optimal, and how to construct the corresponding randomization scheme? We show that: (i) for smooth estimators, antithetic randomization with pairwise correlation $\rho=-1/(K-1)$ is necessary and sufficient for the reducible variance due to randomization to remain bounded as the bias vanishes; (ii) a general construction yields a class of antithetic schemes, within which the jointly normal scheme is minimax optimal; and (iii) for non-smooth estimators with finitely many jump discontinuities, antithetic randomization improves the asymptotic rate of the reducible variance, while a simple control variate restores bounded variance when the discontinuities are known.

Two-way contingency tables are tested for independence throughout applied statistics (genomics, network and text co-occurrence, pharmacovigilance, ecology, survey cross-tabulation), and the routine test reads asymptotic Gaussian tail probabilities off the table cell by cell. On the tables people actually analyze this is badly miscalibrated: across 5,543 real public two-way tables, two-thirds have counts small enough or margins heterogeneous enough that the asymptotic maximum-cell independence scan false-positives at a mean 33% against a 0.05 target, while an exact margin-conditional anchor holds at about 1%. The failure is sharpest when the departure is sparse, the association concentrated in a few cells, where the omnibus chi-square is inefficient and the sharp object is the detection boundary. In three parts we answer where a sparse signal can be seen, which combiner attains that limit, and how to calibrate at small counts. Part I specializes the sparse-detection theory of Ingster (1997), Donoho and Jin (2004), and Chhor, Mukherjee and Sen (2024) to independence: the table is a heteroscedastic Gaussian sequence whose variance profile is the expected-count table, fixed by the margins via CVe^2 = (1 + CVr^2)(1 + CVc^2) - 1, giving an exact signal map and a closed-form separation radius valid under a moderate-deviation growing-count condition. Part II shows higher criticism, with its closed-form Jaeschke-Eicker null, attains that boundary adaptively over unknown sparsity. Part III measures the small-count calibration failure (size 0.3 to 0.6 even under uniform margins, so the driver is the per-cell tail), removes it with exact margin-conditional laws (size 0.003 to 0.009), and gives a routing rule sending large-count tables to the asymptotic gate and small-count tables to the exact anchor. Every claim is reproduced from openly deposited, deterministically seeded code.

We study the problem of testing independence between the coordinate processes of a jointly stationary ergodic binary process from finite observations. We prove that no test is pointwise consistent in probability: any procedure whose power tends to one against every dependent law must have nonvanishing false-positive probability on some independent stationary ergodic law along infinitely many sample sizes. Quantitatively, for some jointly stationary ergodic law with independent coordinate processes, the false-positive probability has limsup at least $1/2$. Thus, even under stationarity and ergodicity, independence cannot be consistently decided from increasingly long finite samples. The proof combines a diagonal construction with rare markers that create detectable dependence at selected scales while converging to a limiting process whose coordinates are independent.