Papers

In the tensor principal components analysis (tensor PCA) problem, we are asked to recover an unknown vector in $\{\pm1\}^N$ from a noisy observation of its $k$-th tensor power, for a fixed $k\ge3$. Tensor PCA exhibits a well-studied information-computation gap: global algorithms such as spectral methods and sum-of-squares work above an algorithmic threshold, below which we know matching lower bounds (against these and some other restricted models of computation). In contrast, local algorithms such as Markov chains were only known to succeed from a random initialization at a polynomially larger signal-to-noise ratio, and were believed to be inferior to global algorithms for tensor PCA. We show that Glauber dynamics, when supplemented with annealing and ergodic averaging, achieves polynomial-time strong recovery above the algorithmic threshold from a uniformly random initialization, demonstrating that local Markov chains can succeed all the way down to the algorithmic threshold from an uninformative cold start.

Resampling is the standard remedy against weight degeneracy in particle filters, yet its effect on the spectrum of the particle second-moment matrix is poorly understood. We study a single resampling step for a d-dimensional ensemble of N particles with d/N -> c. For any unbiased scheme, the resampled empirical second-moment matrix dominates its pre-resampling counterpart in conditional convex order, while the conditional mean is preserved; the excess second spectral moment equals the Frobenius energy of the resampling perturbation. For multinomial resampling, the empirical count law converges to the mixed-Poisson law determined by the limiting weight law. For residual resampling, the count law is determined by the joint limit of floor-fraction pairs: we classify all residual limits of weight sequences converging to delta_1 in law as a one-parameter family between the deterministic identity and the full multinomial count law. Each count law feeds a common covariance-side spectral map, yielding exact consequences: a closed-form zero-eigenvalue mass, a rank-and-spectral-gap transition with a universal lower edge benchmark, computable regular edges, unbounded right support, and scheme-specific second-moment costs. The effective sample size determines one bulk statistic but cannot determine the limiting spectrum. Pre-registered numerical experiments confirm the theory at every level.

Learning interaction structures from spatiotemporal event data is a key statistical problem in fields such as neuroscience, epidemiology, and seismology. We study this problem for Hawkes processes on a continuous spatial domain, where events may trigger future events at different locations. Following the framework developed in “Spatiotemporal Hawkes processes with a graphon-induced connectivity structure” by Baars, Laeven, and Mandjes (2026), we model the spatial interaction mechanism using an underlying graphon that characterizes how the intensities' excitation strength varies across locations. By expanding the baseline intensity, temporal excitation kernel and graphon in orthogonal bases, we represent the problem as a high-dimensional parametric inference problem. We propose a sparse penalized least-squares estimator that recovers the effective interaction structure while controlling model complexity. Our analysis establishes non-asymptotic error bounds in a high-dimensional regime where the number of parameters may grow with the observation horizon, without requiring the conditional intensity to be almost surely bounded. We further establish sparsistency results, showing that the latent interaction network can be consistently recovered. Simulation studies demonstrate the good finite sample performance of the proposed method and an application to large-scale neuronal spike data, comprising 1.81 million events, showcases its practical effectiveness.

We study the recovery of the hidden vertex correspondence between two correlated dense random geometric graphs. Each graph is generated from $n$ latent unit vectors in $\mathbb R^d$; the vectors $u_i^A$ and $u_i^B$ of corresponding vertices satisfy $\langle u_i^A,u_i^B\rangle=\rho$, and the vertex labels of the second graph are permuted by an unknown permutation $\pi$. For every fixed correlation $\rho\in(0,1]$ and every edge density bounded away from zero and one, we give a randomized polynomial time algorithm that recovers $\pi$ exactly with probability tending to one whenever $d\ge(\log n)^L$, where $L$ depends only on the correlation and the density bounds. The algorithm assigns to each vertex a signature formed by counting decorated paths with a constant number of seeds. The number of such seeded templates grows exponentially with the path length, with an arbitrarily large constant base, so the resulting statistics amplify an arbitrarily small constant correlation. We control the dependencies created by the latent geometry in two ways: when $d$ is larger than $n/(\log n)^{L_0}$ for a suitable constant $L_0$, by comparing local moments with those of correlated Erdős–Rényi graphs; for smaller $d$, by spectral estimation of the latent geometry combined with a harmonic correction. We also provide low-degree evidence for a tradeoff between the correlation and the computational cost.

The multivariate linear mixed model combined with restricted maximum likelihood estimation of (co)variance components and solving the mixed model equations is widely used for the analysis of plant breeding data. Several algorithms for restricted maximum likelihood estimation have been proposed, with the average information algorithm being widely used for its computational efficiency, quadratic rate of convergence, and ability to handle large amounts of parameters. A key reason for the scalability of linear mixed models is the exploitation of sparse matrix techniques. Recently, decreasing costs of high-throughput data collection have led to high-dimensional sets of features characterizing genotypes and environments. A convenient way of using these features is through the use of kernel-based covariance structures. However, while convenient, these kernels typically have dense inverses meaning optimizations based on sparse matrix techniques are not always applicable. For these settings, the introduction of general-purpose computing on graphics processing units brought significant computational advantages for dense linear algebra operations. We present a hybrid GPU/CPU implementation of the average information restricted maximum likelihood algorithm. We eliminate sparse equations using supernodal Cholesky factorization and rely on CUDA GPU acceleration to solve the resulting reduced dense system. We show the significant computational advantages of this approach for fully dense problems, as well as problems with a mix of dense and sparse equations.

Bayesian inference for latent point process intensities often relies on per-dataset sampling or optimization, making repeated inference computationally expensive. We propose a neural-network-based amortized variational framework that learns a reusable inference rule mapping a variable number of replicated point patterns to a distribution over their shared latent intensity. Trained from event times using the point process likelihood, the model learns an amortized inference rule without ground-truth intensity supervision. At test time, it performs inference of the latent distribution in a forward pass. On synthetic data, additional realizations improve intensity recovery and yield narrower credible intervals with broadly stable coverage. Our method is tested on bicycle-sharing data and achieves the highest mean held-out log likelihood among the evaluated methods while requiring substantially less inference time than the Log-Gaussian Cox Process baseline.

Non-exchangeability and radial asymmetry have got substantial attention in the literature, but relation between the two has not been studied yet. We propose the question: if a copula has a given degree of radial asymmetry, what is the range of possible values of the measure of non-exhangeability. We determine the set of all points whose first coordinate is the radial asymmetry of a copula and the second coordinate is its non-exchangeability. We show that this set is a non-convex hexagon.

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.

Modern deep classifiers trained by gradient descent need many passes over the data to become accurate, and even then are not always well calibrated; post-hoc temperature scaling corrects calibration after training but cannot influence training-time accuracy or convergence. We propose a variational filtering algorithm for the last layer of deep classifiers, whose parameters are treated as the latent state of a linear-Gaussian state-space model. Each mini-batch is assimilated by an exact prediction step followed by a closed-form variational update, made conjugate by an independent binary probit (IBProbit) surrogate for the categorical likelihood. The filter's process noise stabilises the updates under feature drift and keeps the predictive distribution from collapsing to a point estimate, allowing the same procedure to be used for linear probing on a frozen backbone and for full-network fine-tuning. We evaluate on DINOv3 backbones against two baselines tuned under an identical protocol: AdamW with cross-entropy and post-hoc temperature scaling, and Variational Bayesian Last Layers (VBLL). As a linear probe on frozen features, IBProbit attains the best accuracy and negative log-likelihood on all five downstream tasks at the largest backbone, and on two of the five at the smallest. In out-of-distribution detection with CIFAR probes on the largest backbone, IBProbit reduces the false-positive rate at 95% true-positive rate by a factor of 1.4 to 3.2 relative to the better of the two baselines, on all three dataset pairs. When the whole network is fine-tuned on ImageNet-1k, IBProbit reaches the final performance of the cross-entropy baseline in half the wall-clock time.

Accurately assessing the convergence of Monte Carlo Markov Chains is crucial for ensuring their validity in Bayesian inference problems. However, convergence assessment in high dimensions is known to be difficult, due to most convergence diagnostics becoming insensitive beyond a few dozen parameters. Conventional convergence statistics analyze marginal distributions which in high dimensions neglects increasing levels of information from the joint and conditional probability distributions. Thus we propose $\texttt{NFDiag}$, which taps into the information of the full multivariate joint posterior with a normalizing flow (NF). Our convergence diagnostic thereby also analyzes information from conditional distributions. In other words, by working in the full multivariate space, it is sensitive with respect to the geometry of the chain in high dimensions. Our diagnostic utilizes that NFs assume independent and identically distributed (iid) training samples, which is a testable assumption that will be violated by a non-converged or strongly correlated chain. Our diagnostic trains NFs on separate parts of a chain, and then assesses whether these correspond to the same multivariate target distribution. We test our diagnostic in up to 250 dimensions, and demonstrate it is universal and unbiased, and detects problematic chain quality also when Gelman-Rubin $\hat{R}$ and the multivarite effective sample size do not.

In a variety of multi-agent and physical environments, an element's update depends strongly on a few other elements and weakly on many, and both the identity and the number of the elements that matter change with the state of the system. Such graded sparse mechanisms are common, for instance, in flocks, crowds, learned particle simulators and the neighborhood aggregation of graph networks, yet they have received little systematic study as a class of functions with the grading given. In this paper we study the supervised learning of the mechanism when the grading is provided by the environment, as a weight on every element. First, we define the class by axioms, where each element's influence is capped by its weight, the dependence on the weights is regular, and no element is distinguished, and prove a representation theorem: the class is exactly the Lipschitz ball of a transport metric on weighted multisets, and the metric is forced by the axioms. Next, the geometry of the representation gives generalization bounds that depend on the data in two ways. One is via the sparsity profile of the data – the law of how many elements matter and by how much – which bounds the rate over all designs with that profile, with an exponent from the profile and no dependence on the number of elements. In a number of natural cases we also show matching lower bounds, up to constants in the exponent. The other is empirical: due to the representation result, classical empirical discrepancy bounds apply to the class as they stand, and the resulting bound is the transport distance between the two halves of the sample in the class's own metric, computed exactly by a matching, with an additional concentration term. Finally, we study a simulated Vicsek flock under four interaction rules at two noise levels, and show that the certificate is non-vacuous, a small fraction of the trivial value on every rule.

We study the last iterate of stochastic gradient descent (SGD) on streaming least squares with discretely distributed features (Slock's model), which decouples the dynamics across eigendirections. Under power-law source and capacity conditions, we derive the exact finite-time risk and asymptotic equivalents of constant and Warmup-Stable-Decay (WSD) schedules. With a tuned base learning-rate, the resulting scaling law decays as $\Theta(T^{-q})$, where $T$ is the horizon, and $q$ is independent of the cooldown length $c \in (0,1]$ but is strictly smaller for the constant schedule: any cooldown of fixed relative length $c > 0$, however short, improves the risk by a polynomial factor, while its length only affects the constant, which is minimized by linear decay ($c = 1$). Finally, we show that the optimal risk depends on the batch size $b$ and horizon $T$ only through $bT$, up to a critical batch size.

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.

The estimation of heterogeneous treatment effects is a central objective of causal inference. However, state-of-the-art algorithms to solve this task are based on machine learning principles, and thus behave as black boxes. Consequently, the identification of covariates driving heterogeneity remains a challenging and open problem. To address this limitation, we introduce CausalTreeHFD, an interpretable extension of causal forests, based on the Hoeffding functional decomposition adapted to tree ensembles. Our approach decomposes predictions of conditional average treatment effects (CATE) into contributions from individual covariates and pairwise interactions, providing a direct characterization of the factors driving treatment effect heterogeneity. We establish the consistency of CausalTreeHFD, which converges towards the Hoeffding decomposition of the underlying CATE function, under standard assumptions for the analysis of causal forests. Experiments on simulated and real datasets show that CausalTreeHFD accurately identifies relevant covariates and interactions, and precisely recovers their contributions to heterogeneous treatment effects.

We study Fisher information about the amplitude of a known smooth mean shift from the first Dirichlet eigenvalue and spatial responses of a Liouville random medium. The complete response record determines the speed measure and the Gaussian free field, whereas every finite record strictly loses information; an optimal design exists in each finite-dimensional probe space. We quantify this gap as the known Liouville parameter $\gamma$ tends to zero. If the probe space contains an explicit matching probe, the minimum loss equals $\gamma^2$ times a positive compression error for a deterministic second-chaos operator, up to $o(\gamma^2)$, with designs allowed to depend on $\gamma$. After enlarging the probe dictionary, the optimal $m$-response coefficient is asymptotic to $\|Lg\|_{L^1(D)}^2/(8m)$, where $L=-(2\pi)^{-1}\Delta_D$ on the unit disk $D$. A quadratic correction of the readouts attains the same fixed-design loss to second order.

Realised kurtosis measures intraday tail risk but compresses the shape of the return distribution into a single scalar, so trading days with similar kurtosis may differ in how probability mass is spread across the peak, the shoulders and the tails. We propose a geometric description of intraday return shape based on the modified Brownian motion model (MBMM), which deforms a Gaussian benchmark through a nonlinear feedback governed by a shape parameter and a distortion intensity. Our main contribution is a peak-shoulder-tail decomposition that measures where the fitted return density gains or loses probability mass relative to the Gaussian, complemented by two feedback landmarks that explain how the shape parameter generates these deformations. We estimate the model at the stock-day level using a bin-free simulation-based minimum-distance estimator that matches empirical and model-implied intraday quantiles, and we document its finite-sample reliability by Monte Carlo. Using regular session one-minute returns for ten liquid U.S. stocks from May 2021 to April 2026, we find that the MBMM provides a substantially better in-sample fit than geometric Brownian motion on every estimator-ready stock-day, and that days with closely matched realised kurtosis can display materially different peak-shoulder-tail geometries. Finally, we show that the estimated shape parameters behave as persistent daily state variables: they are strongly forecastable relative to a random-walk benchmark, though their predictability reflects persistence around a slowly moving level rather than rich autoregressive dynamics. The proposed decomposition offers a risk-shape measurement layer that complements realised kurtosis by describing not only how large intraday tail risk is, but how it is formed, giving investors and risk managers additional distributional information for risk monitoring and, potentially, prediction.

We propose an anomaly detection method for high-dimensional mixed-type data that maps heterogeneous features to a common scale via standardized moment-based transformations, yielding interpretable feature-level anomaly scores. These scores are combined through a dependence-aware aggregation step to construct the Mixed-data Influence (MDI) score, accounting for correlations across features and providing a unified measure of how strongly an observation departs from the reference data distribution. By incorporating higher-order moments, the MDI can detect anomalies that arise not only from location shifts, but also from changes in variability, skewness, and other distributional characteristics. The aggregation mechanism allows our method to capture both sparse anomalies concentrated in a small number of features and weak, coordinated deviations distributed across many features. We empirically investigate the distribution and thresholding behavior of the MDI score across different sample sizes, dimensions, dependence structures, and data types. Simulation experiments and two real-data applications, to multi-omics and ESG data, evaluate finite-sample calibration, detection performance, and feature-level interpretability.

We study sampling from $\pi(x)\propto\mathrm{e}^{-V(x)}$ on $\mathbb{R}^d$, where $V$ is $\alpha$-strongly convex and $\beta$-smooth with condition number $\kappa=\beta/\alpha$. For this problem, the proximal sampler preserves the target exactly and recovers overdamped Langevin dynamics in the small-step limit, but its Gaussian resampling step discards directional information. Underdamped Langevin dynamics retains such information through momentum to accelerate exploration, motivating a proximal construction that combines target preservation with directional persistence. We introduce reflected proximal dynamics, which combines proximal sampling with an exact reflection of the auxiliary variable, conditional kinetic motion, and controlled refreshment. We develop two realizations: a flux-resampled underdamped half-turn and a persistent-momentum quarter-turn. From a warm start with constant $\chi^2$ divergence, both methods achieve total variation error $\varepsilon$ using $\widetilde O(\sqrt{\kappa}d^{1/4})$ expected oracle queries. We further remove the supplied warm-start assumption by combining either sampler with existing posterior samplers in a proximal initialization stage with an optimized coarse step size. Accounting for both initialization and sampling yields $\widetilde O(\kappa^{6/7}d^{1/4})$ expected oracle queries. In both regimes, the dependence on the target accuracy $\varepsilon$ is $\mathrm{polylog}(1/\varepsilon)$, making reflected proximal dynamics a promising framework for accelerated high-accuracy log-concave sampling.

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.

A compact metric tree is the natural domain of spatial data recorded along a river network, a root system, or the branched processes of a cell, and a suitable distance between two points is the length of the path joining them. That distance admits almost none of the covariance models of spatial statistics, since only the roughest members of each family stay nonnegative definite, and the resulting fields are nowhere mean square differentiable. The repair proposed here is a change of metric. Comparing two points through their distances to the leaves gives a Hausdorff-type metric built from the shortest-path one. The new metric has a closed form, agrees with the shortest-path one inside every edge, and is strictly smaller across branching points. Any correlation function valid in a Euclidean space of dimension equal to the number of leaves may then be composed with it, returning the Matérn, powered exponential, and generalized Cauchy families at every degree of smoothness, and, on trees with few leaves, also correlations that change sign or have compact support. The same construction yields the multivariate Matérn and space-time models of the Gneiting class. In simulations covering a range of smoothness and of spatial scale, maximum likelihood provides estimates without appreciable bias, and on a reconstructed microglial cell the smooth models reduce the cross-validated prediction error by half.

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.

Existing cumulative regret guarantees for Bayesian optimization (BO) do not necessarily characterize the quality of the algorithm's final output. This gap is critical under noisy evaluations, where the goal is to return a single high-quality recommendation. To address it, we introduce efficient exploration strategies that can be incorporated into any BO procedure by occasionally replacing the acquisition query with an exploration query. To validate this approach, we bound the Gaussian process posterior standard deviation over the search space via the fill distance of the queries, which measures the largest region left unexplored. We then use this result to derive upper bounds on the inference regret, i.e., the gap between the global optimum and the function value at the algorithm's returned point. Empirically, these exploration strategies improve the final recommendations of existing BO algorithms on synthetic benchmarks, hyperparameter optimization, and a control task.

Learning equilibria from data is difficult when several decision makers interact indefinitely and the system state has no prescribed bounds. We examine discounted general-sum stochastic games with compact action sets and an unknown nonlinear drift observed through Gaussian transition noise. The proposed procedure uses a single uninterrupted trajectory to update Gaussian-process models of the drift and to solve successively refined finite games. We prove a perturbation result showing that stability of equilibrium selection, together with uniform KL accuracy of Gaussian transition models, is sufficient for regional uniform convergence of the induced values and approximate Nash behavior in the original game. The region specifies where the guarantees are evaluated; strategies and dynamics continue to act on the entire Euclidean state space. For GP estimators, simultaneous confidence bounds yield a rate controlled by posterior uncertainty, kernel information gain, and joint-state dimension. Numerical studies illustrate how kernel choice and dimension appear in the predicted convergence behavior.

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.

The $r$ largest observations of a block form a dependent random vector whose support moves with all three parameters of the $r$-largest generalized extreme-value ($r$GEV) model. For a fixed number $r$ of order statistics, we show that the model is differentiable in quadratic mean if and only if $k_0<1/2$, where $k=-\gamma$ is the shape parameter; the threshold is thus the same as for block maxima, because only the largest or the smallest retained order statistic can reach the moving boundary. In this range the Fisher information is positive definite and decomposes into the block-maximum information plus contributions of the lower order statistics, and it equals the negative expected Hessian, which we prove without using any estimator. The model is then locally asymptotically normal, whereas for $k_0\ge1/2$ root-$n$ local asymptotic normality fails. Numerically, using the five largest order statistics instead of the block maximum alone lowers the efficient asymptotic standard error by about $45\%$ for the shape parameter and by $20$–$24\%$ for the 100-year return level in heavy-tailed cases ($k_0\le-0.2$), but by only about $2\%$ for the return level when $k_0=0.4$.

Real-world public health policies are often enacted by local jurisdictions, but their effects can ripple through connected populations. For causal evaluation, this creates a difficult spatiotemporal problem: intermittent policy updates create path-level treatments, and because spillover depends on neighboring treatment paths along with the focal-unit's own history, controlling time-varying confounding requires high-dimensional, multimodal information from both sources. To address this challenge, we develop a neural spatiotemporal spillover framework that extends potential-outcome causal analysis to networked policy-systems with unit-level effects and structured interference. The framework indexes each unit's potential outcomes through own-treatment and neighborhood spillover-exposure paths, yielding direct, spillover, and total effects. We introduce an LSTM-based representation-learning encoder that integrates each unit's own and neighborhood histories into compact sequential summaries for downstream modeling and causal estimation. These summaries feed a self-normalized kernel inverse-probability-weighting estimator for continuous spillover exposure. Beyond theoretical guarantees, simulations show improved performance over outcome-regression baselines; in a practically motivated county-level COVID-19 policy study, the framework finds both business-and-economic restrictions and education-and-childcare policies associated with lower subsequent case burden, with larger estimates for business-and-economic restrictions. Overall, it provides a scalable, policy-interpretable framework for evaluating unit-level effects in high-dimensional, multimodal spatiotemporal systems when interventions spill over across space and time.

A quickest change detection problem is studied where the decision-maker has access to multiple experiments, but can choose only one experiment at a time to detect the change. The true post-change distribution under each experiment is unknown, but known to belong to a family of distributions. It is assumed that each of these families of distributions has an element that is least favorable in a well-defined sense. An algorithm is proposed that is designed using the least favorable distributions. It is shown that the algorithm is first-order asymptotically optimal, as the rate of false alarms goes to zero. This algorithm minimizes the worst-case average detection delay over all post-change distributions within the assumed families, subject to a constraint on the average run length to false alarm and constraints on the average fraction of time each experiment is performed before change.

Backward conformal prediction (BCP) fixes the size of a prediction set and estimates its miscoverage, but assumes calibration and test data are exchangeable. We extend BCP to covariate shift through importance weighting, preserving its guarantees when the density ratio is known. We show that BCP's tight miscoverage report depends only on the relative weight assigned to missed calibration examples, explaining why a better domain classifier can sometimes produce a worse estimate. Motivated by the fact that BCP depends only on nonconformity scores, we propose estimating density ratios on sorted score vectors rather than high-dimensional inputs. We prove that this approach recovers the same report asymptotically under an interpretable condition and characterize the resulting bias otherwise. Experiments with simulated and benchmark distribution shifts show that the proposed method improves BCP's miscoverage estimates, while score-based weighting consistently yields more accurate estimates than embedding-based weighting.