Papers

In this article, we explore a new paradigm for statistical inference. The approach centers around the point estimate based on observed data, simulating replicates using the estimate as the truth to produce clones of the estimate, with inference deriving from the clone distribution. It avoids prospective finite-dimensional model assumptions, but it makes no probabilistic claims concerning the truth; it suggests an alternative system of uncertainty quantification that is operable in nonparametric function estimation. The procedures are demonstrated using examples of smoothing spline ANOVA models in nonparametric regression. The paradigm also applies in parametric regression, where the proposed inference closely resembles traditional inference operation-wise. Conceptual discussions are scattered throughout.

Many questions across the sciences take the same form: several coupled series are observed together, and the analyst wants to know not merely that they move together but which one moves first, and how strongly. This paper sets out a complete method built on one organising idea: the direction of a coupled system is exactly the part of its behaviour that changes when the record is played backwards. Tools built on contemporaneous covariance alone (correlation matrices, distance measures, spanning trees, undirected centralities, principal components) carry no information about direction: a reversible system and a circulating one can share identical covariance at every sampling of the same point-in-time record. Formally, direction is a circulation matrix carried by the lagged covariance. Its vanishing is exactly statistical time reversibility for linear systems, feature maps carry the characterisation to nonlinear ones, and under the Gaussian benchmark its magnitude is an entropy-production functional of the identified circulation, the quadratic component of the divergence per unit time between the forward and reversed records. Around this estimand we build a cross-fitted estimator removing first-order bias, delete-block jackknife standard errors, and a randomisation test exact under its stated block null, with a familywise correction and a nonlinear extension. A sampling theory says when the arrow is measurable at all, and a design layer separates transmission from the ordering of clocks. A laboratory of four systems with known answers compares the method with correlation networks, Granger causality, transfer entropy, and connectedness indices, reporting the failures of each when read as a measure of direction, including our own. Complete algorithms and worked examples in two languages make the paper the base reference for a series of applications.

Early-phase clinical trials for dose selection typically enrol few patients and aim to identify doses that are both safe and promising for further study. While traditional approaches identify the maximum tolerated dose, modern trials for targeted therapies often seek the optimal biological dose, defined as the lowest dose achieving sufficient biological activity with acceptable safety. In immunology settings, assessment of biological activity is based on multiple biomarkers or clinical endpoints. Clinicians leading dose-selection efforts would thus benefit from transparent summaries of the probabilities of observing combinations of biomarker outcomes across doses. However, such inference is challenging in small samples where complex modelling assumptions are difficult to verify. To address this limitation, we propose COBRA-DOSE, a framework for posterior predictive inference based on two endpoints that models dependence via copulas and accounts for uncertainty in both marginal distributions and dependence structures through Bayesian model averaging. This approach avoids reliance on a single model and yields interpretable quantities for clinical decision making. We demonstrate the performance of COBRA-DOSE using DEN-181, a phase I immunology trial in rheumatoid arthritis. We also provide a general implementation of our approach through the CobraDose package in R.

Given a dataset where a portion of the samples are contaminated, our goal is to recover the underlying clean population distribution. To this end, we propose Wasserstein Filtering (WF), a novel sample selection framework that discards a fraction of suspicious samples and estimates the target distribution using the empirical measure of the remaining data. The core insight is to select a subset of samples whose empirical distribution maximizes its Wasserstein distance to the fully contaminated empirical distribution, thereby preferentially isolating and removing geometrically influential outliers. To render this optimization computationally tractable, we introduce three algorithms: a marginal screening scheme, SinkMarg, and two joint optimization algorithms, SinkWF and SlicedWF, leveraging entropic optimal transport and sliced Wasserstein approximations, respectively. On the theoretical front, we introduce the Far Exclusion and Local Projection (FELP) contamination model, which characterizes corruptions consisting of well-separated outliers and locally indistinguishable perturbations. Under this model, we prove that the WF estimator achieves minimax optimality over distribution families with bounded covariance. Extensive numerical experiments on synthetic datasets, benchmark anomaly detection suites, and robust generative learning with diffusion models demonstrate that WF serves as a highly practical, model-agnostic preprocessing tool. It delivers competitive outlier detection performance and provides substantial downstream benefits for generative modeling under heavy contamination.

Cardiac fibrosis reduces electrical conductivity and is a leading cause of arrhythmia. Arrhythmic waves typically rotate around non-conducting fibrotic patches, so the geometry and topology of these patches (spatially isolated regions of fibrotic tissue within the heart muscle) play an important role in arrhythmia dynamics. Despite their clinical relevance, these structures remain poorly understood. We address this open problem using histopathological images of human hearts affected by cardiac fibrosis. Each patch is represented as a spatial graph via skeletonization, where nodes are embedded as points in Euclidean space and edges encode geometric properties of the underlying tissue. The core methodological contribution of this work is the introduction of spatial graph space, a metric space equipped with a rotation-invariant Fused Gromov Wasserstein metric that enables comparison of spatial graphs with differing numbers of nodes and edges. Building on this, we perform a distribution-level statistical testing and depth measures for spatial graphs. To enable the interpretation of the spatial graph sample distribution, we introduce DepthPlot, a novel visualization tool for depth measures in metric spaces. Applying our methodology to compare patients and the spatial position of patches within the ventricles, we find that fibrotic textures exhibit strong patient-specific features, while some hearts display notable geometric similarities, potentially reflecting shared pathological mutations or other unknown factors. Through quantitative depth measures, we characterize test outcomes via central and peripheral spatial graphs, demonstrating that the proposed framework yields statistically and clinically meaningful insights into fibrotic texture characterization.

The construction and the stability of bridges between mixture probability measures are investigated. Given a collection of Markov transitions bridging the components of two mixtures, together with a coupling of the mixture labels, we construct a bridge between the mixtures on an extended state space, and we show that this bridge inherits the entropic optimality of the component bridges. We then turn to Gaussian mixtures bridged by Gaussian Schrödinger bridges, whose parameters are generally unknown and must be estimated. We prove a $2$-Wasserstein continuity theorem between the exact bridge and its plug-in approximation. The analysis rests on a representation of the Riccati fixed point map in inverse coordinates; this map is globally $1$-Lipschitz and requires no matrix inversion. These results are applied to three estimation schemes, with $d$ the dimension and $N$ the number of samples. For Gaussian mixtures estimated by the Expectation-Maximization algorithm, the squared $2$-Wasserstein error is of order $(d\log N/N)^{1/2}+d^2\log N/N$ with probability at least $1-10N^{-1}$. For single Gaussian marginals estimated by their sample moments, it is of order $d^2/N$ in expectation. For the regularized empirical Monge map based on an $\epsilon$-inflation of the sample covariance, the mean squared error is of order $d^2\{(1+\epsilon^{-1})/N+\epsilon^2\}$; as soon as $N$ is a sufficiently large multiple of $d$, it reduces to $d^2/N+d\,\epsilon^2$, up to an exponentially small term. These estimates are illustrated numerically.

We prove new comparison results between the Wasserstein distance and its sliced and max-sliced counterparts. First, we show that the Hölder exponent~$\frac{2}{d+2}$ obtained by Bobkov and Götze for the max-sliced 1-Wasserstein distance on the unit ball is optimal for every $d \geq 2$, settling a question raised in their work. Second, we show that sharper comparisons are possible under stronger structural assumptions: if $\nu$ is a discrete measure and the optimal coupling between $\mu$ and $\nu$ transports each point to a nearest atom of $\nu$, then $W_p(\mu, \nu) \leq C \sqrt{d}\, K \, \mathrm{SW}_{p,1}(\mu, \nu)$ for a universal constant $C$, where the complexity parameter $K$ is always at most the number of atoms $N$ and can be substantially smaller. This complements a similar bound due to Park and Slepčev. An analogous bound holds for the sliced Wasserstein distance based on $k$-dimensional projections. Finally, using a construction from geometric discrepancy theory due to Chen and Travaglini, we prove that the linear dependence on $K$ in this bound cannot be improved, up to polylogarithmic factors.

This study develops a weighted framework for measuring the discrepancy between two nonnegative lifetime distributions through cumulative past extropy. We propose two measures, referred to as the weighted cumulative past extropy inaccuracy (WCPEI) and the weighted cumulative past extropy Kullback-Leibler divergence (WCPED). The generalized weight function is considered in this study. We investigate a number of theoretical properties of these measures. Empirical distribution function-based nonparametric estimator is subsequently constructed for the weighted cumulative past extropy inaccuracy ration (WCPEIR). Its finite-sample behavior is studied through Monte Carlo simulation experiments for different sample sizes. To illustrate the practical relevance of the WCPED, two applications are considered. First, an extropy-based goodness-of-fit procedure for testing uniformity is developed using the proposed divergence measure. Its power is then compared with that of several established uniformity tests under a variety of alternatives. Second, an image analysis application is presented in which the proposed measure is employed to assess changes in the distributions of pixel intensities when the image resolution is altered. The framework is further extended to a dynamic setting by conditioning on the lifetime information available up to a specified time point. This leads to the dynamic weighted cumulative past extropy inaccuracy (DWCPEI) and dynamic weighted cumulative past extropy divergence (DWCPED). Their theoretical properties are derived, and the corresponding nonparametric estimation procedures are proposed. The finite-sample performance of these estimators is evaluated through simulation studies using R software.

We develop bootstrap-assisted robust binary segmentation (BARBS), a recursive binary segmentation method for multiple change point detection under general nonstationary temporal dynamics. A novel Gaussian multiplier bootstrap for the CUSUM statistics is proposed, offering robustness to complex dependence structures. Through meticulous calibration of the critical values at each stage of the recursion, BARBS ensures control of the Type I error under the null hypothesis of no change points. When change points are present, BARBS identifies the correct number of changes with a prespecified probability, and the resulting change point location estimators attain the same uniform consistency rate as classical binary segmentation. Building on this, we introduce second-stage refined estimators that achieve the optimal individual localization rate, and establish their asymptotic distributions and nearly optimal uniform localization rates under both fixed and vanishing jump magnitudes. Extensive numerical experiments across various settings confirm the robustness and superior performance of BARBS relative to existing approaches. To illustrate the practical relevance of the proposed methodology, we analyze U.S. inflation data, yielding change points that align with several documented macroeconomic episodes.

We propose a distributed selective inference framework tailored for high-dimensional quantile regression. To enable valid post-selection inference in this context, we address the computational challenge posed by the non-smooth quantile loss via a response-surrogation strategy. This strategy transforms the problem into a penalized least-squares formulation, thereby facilitating distributed selective inference. For valid post-selection inference, a randomized procedure is introduced, in which the Lasso selection event is characterized through the associated Karush-Kuhn-Tucker conditions and the conditional distribution of the aggregated estimator is derived given the selection event. The resulting algorithm requires only three rounds of communication between local machines and the central server. Under standard regularity conditions, we establish the asymptotic validity of the proposed procedure and develop a large-deviation approximation to the selective likelihood for computationally tractable implementation. Simulation studies and a real-data application demonstrate the satisfactory finite-sample performance of the proposed method.

Estimating the distribution of relaxation times (DRT) fromelectrochemical impedance spectroscopy (EIS) is an ill-posed inverse problem that is highly sensitive to regularisation choices. We propose a physics-informed convolutional autoencoder that estimates DRT directly from EIS data without spectrum-specific tuning. A discretised relation between impedance and the DRT is embedded in the training process, constraining the network to produce impedance-consistent distributions. The model resolves overlapping relaxation processes in synthetic two-ZARC spectra and accurately reconstructs measurements from three independent solid oxide fuel and electrolysis cell datasets, with range-normalised errors below 1.1%. Decoder-probe analysis shows that the learned latent representation is organised according to relaxation timescale. Distances in this latent space capture operating changes, hydrogen-shortage events, and long-term degradation. The same lightweight architecture is applied across all datasets without modification, providing consistent DRT estimation and an interpretable basis for condition monitoring.

Spherically invariant (SI) random processes can model impulsive noise and unreliable measurements. Recently, the mixture noise of Gaussian and SI components has been used in deterministic maximum likelihood direction finding. In this context, the Expectation-Conditional Maximization (ECM) algorithm, an extension of the expectation-maximization algorithm, has been applied and designed. However, simulation results show that the ECM algorithm always improperly converges. In this article, the ECM Either (ECME) algorithm, an extension of the ECM algorithm, is applied and designed, which additionally utilizes the actual log-likelihood function to first update partial parameter estimates at every iteration and does not need to initialize all parameter estimates. Moreover, the deterministic Cramer-Rao low bounds (CRLBs) of DOA estimators are derived and compared. Simulation results indicate that the ECME algorithm exhibits proper convergence and its root mean square errors of DOA estimates asymptotically approach the CRLBs as the signal powers increase, i.e., the derived CRLBs are correct.

Causal mediation analysis is typically formulated under no interference, an assumption often violated in networked populations. We develop a nonparametric framework for a single large observed network that allows simultaneous treatment and mediator spillovers and high-dimensional network confounding. Exposure and mediator mappings define causal estimands without restricting the true interference mechanism, separating own from spillover effects without prespecified aggregation models. Under strengthened conditional independence conditions, we identify own controlled direct, natural direct, and natural indirect effects and give primitive sufficient conditions in terms of structural errors. We construct augmented inverse probability weighted estimators that are doubly robust for controlled effects and multiply robust for natural effects, using graph neural networks to learn high-dimensional nuisance functions from node features and the adjacency matrix. Under approximate neighborhood interference, weak network dependence, and suitable first-stage rates, we establish asymptotic normality of the effect estimators and consistency of a network HAC variance estimator. In simulations the graph neural network estimator outperforms machine learning methods built on hand-constructed neighborhood features, and a reanalysis of an agricultural insurance experiment in rural China finds insurance knowledge to be a substantive mediating channel while perception-based mediators are not.

We present the mathematical foundations of linear independent component analysis (ICA) models based on standard literature in a self-contained note. It is aimed at readers with a background in measure-theoretic probability theory. We first develop the theory of the characteristic functions of probability measures on $\mathbb{R}^d$, including their analyticity and the way in which they determine and characterise the distributions. We then focus on several identifiability results of ICA models with successively strengthened assumptions on the sources: from merely non-constant, to non-Gaussian, to Gaussian-free independent sources. Under the strictest assumptions, we show that the independent sources are identifiable up to translation, permutation, scales and signs, and this even in the presence of additive Gaussian noise. Furthermore, we present the online equivariant gradient descent ICA algorithm for recovering the independent sources from data, in the standard complete noiseless non-Gaussian ICA setting.

Many operational decisions are sequences of interventions under a cumulative resource limit, such as a maintenance schedule within a crew-hour budget. Choosing among them calls for the outcome and the cumulative cost each would produce, counterfactual quantities identified from observational data. Two strategies with the same expected cost can exceed the budget at very different rates, so constraining the mean does not bound how often an overrun occurs. Prior two-step architectures, recently extended to continuous doses, constrain the mean cost rather than its tail and allocate at a single decision point. Methods that do bound a cost tail take its distribution from a specified model rather than identifying it from data. We present a predict-then-optimize framework. In the prediction step, any estimator returning an outcome value and a cost distribution supplies what the decision rule consumes, so the predictor is interchangeable. In the optimization step, a chance-constrained selection over a finite candidate set bounds the probability that the cumulative cost exceeds the budget. That tail does not decompose across stages, so each strategy is scored whole. Sweeping the tolerated violation probability traces a safety-utility frontier, and distribution-free finite-sample bounds cover violation and outcome shortfall. Four of five environments, spanning clinical treatment and equipment maintenance, supply exact counterfactual ground truth; the fifth carries real outcomes from a digital-health micro-randomized trial. Across them, the rule holds the budget where a point-estimate rule overruns it, at an outcome cost the frontier makes explicit. All code is available at https://github.com/mfriendly/counterfactual-chance-selection

We develop the statistical and algorithmic theory of inverse optimal transport (IOT) under the feature-parameterized cost C_theta(i,j) = -theta^T phi(i,j). The core technical contribution is the Sinkhorn linearization – the implicit-function sensitivity of the entropic OT plan to the cost – together with its spectral proxy, a formula that is spectrally exact yet geometrically transparent. The restricted Hessian on the tangent space satisfies the spectral sandwich (pi_min/epsilon) I <= H_T^{-1} <= (pi_max/epsilon) I, yielding the single core bound sigma_min >= (pi_min/(a_max epsilon)) sqrt(lambda_min(Sigma)) that drives the entire theory. On this core we establish four theorems and one observation. T1 (identifiability): theta is globally injective on the quotient of the gauge kernel, with dimension bound F <= (K-1)^2. T2 (sparsistency): the l1-penalized estimator recovers the true support under irrepresentability and score concentration, with exponential failure probability. T3 (well-posedness): the feature-moment map M(theta) = Phi^T x_theta is strongly monotone, and the inverse is Lipschitz with constant L <= epsilon ||Phi^T S_a||_op / (pi_min lambda_min(Sigma)). T4 (convergence): local strong convexity with mu >= pi_min^2 lambda_min(Sigma) / epsilon^2 guarantees monotone gradient descent convergence. O5 (misspecification): the estimator converges to the OT-model projection of the truth; the Holder continuity of the projection map is assessed numerically, yielding setting-dependent empirical exponents alpha_eff in (0,1).

Socioeconomic indicators in England exhibit complex spatial patterns that are not well captured by standard approaches based on averages or broad geographic classifications. We propose a method for comparing areas based on their internal spatial structure, using a multiresolution representation derived from the discrete wavelet transform. The method embeds areal data into a regular grid, extracts local windows, and represents each as a set of scale- and direction-specific detail coefficients. A dissimilarity measure, defined over these coefficients and minimised over rotations and reflections, is used to identify contiguous sets of statistical units with similar spatial structure. We apply the approach to England's 2025 Index of Multiple Deprivation at the lower layer super output area level. We show that areas with similar internal structure are often found across regions, levels of urbanisation, and average deprivation, challenging the use of these categories as proxies for local geography. The results provide a framework for identifying comparable places based on how deprivation is distributed within them, with implications for policy evaluation and transferability.

To avoid missing important variables and their connections in networks, more and more variables are included in network analysis. Here we show that in a setting with many more parameters than observations (high-dimensional) it is possible to get a conservative (i.e., low false positive rate) estimate of the neighbourhood for each node (which connections are in the network). A neighbourhood is often estimated with a linear model, and this leads to two interesting cases: (i) If the true model is linear, then neighbourhood selection work reasonably well, and (ii) if the true model is nonlinear, then neighbourhood selection requires a penalty for the high dimensions. Here we show the impact of the ridge parameter on the mean squared error, and how this leads to low test variance and hence to neighbourhoods with large numbers of edges. We connect these insights with results from machine learning, where the so-called double descent (when more parameters are included than observations, the mean squared error goes down a second time) has put the traditional view on model selection upside down. Essentially, for adequate neighbourhood selection in models with a large number of parameters, the volume of the model space needs to be included in the penalty. Most neighbourhood selection methods (e.g., Lasso, AIC, BIC) lead to spurious edges (high false positive rate), but we prove that in the high-dimensional setting, minimum description length leads to correct neighbourhood selection or smaller (low false positive rates) in both cases when either the model is correctly or incorrectly assumed linear

We study the laws of large numbers for the largest eigenvalues among all principal minors of Wishart matrices and deformed GOE matrices. We propose a new method based on identifying the deterministic sets to which the random sets formed by suitably normalized principal minors converge in Hausdorff distance, thereby reducing the original extreme-value problems to finite-dimensional convex optimization problems. We demonstrate the effectiveness of this method in regimes not covered by the existing second-moment arguments in \cite{cai2021asymptotic,hu2023extreme}. For deformed GOE matrices with fixed minor size \(k\), we determine the limit for every diagonal variance \(a>0\) and identify a phase transition at \(a=2\). Above the transition, the limiting constant satisfies an explicit recursion with no close-form expression, and the optimizers exhibit a nested hierarchical structure, thereby resolving the case left open in \cite{cai2021asymptotic}. For Wishart matrices with general sub-Gaussian entries and fixed \(k\), we characterize the limit through an entropy-constrained deterministic convex set. When the entries are standard Gaussian, we solve the resulting optimization problem explicitly and obtain the exact value of the limiting constant.

We consider the classical additive measurement-error model $X=Y+Z$, where the latent random variable $Y$ has unknown distribution $F_Y$ and the error $Z$ has a known distribution. We develop direct estimators for three functionals of $F_Y$: (i) $F_Y(x)$ at continuity points; (ii) interval probabilities $F_Y(y)-F_Y(x)$ when $x<y$ are continuity points; and (iii) the size of a jump at a prespecified discontinuity. We derive non-asymptotic bias and variance bounds, and establish asymptotic unbiasedness and consistency. Unlike previous work, we do not require $F_Y$ to admit a density, have a mixture representation, or satisfy global Sobolev smoothness assumptions. The framework accommodates arbitrary latent distributions, including those with both discrete and continuous components, and distributions with multiple jumps. These results rely on a link between Fourier inversion theorems and the algebraic structure of a class of estimators proposed in Mynbaev, Martins-Filho and Henderson (2022). A simulation study evaluates feasible tuning procedures and, where available, compares the finite-sample performance of the proposed estimators with existing methods.

In the linear regression model, we construct a nonparametric estimate of the regression parameter vector $\boldgreek\beta$ that is insensitive to a possible nuisance autoregression in the model errors. The main tool for estimating $\boldgreek\beta$ is based on the autoregression rank scores of the model. The resulting estimator is invariant to the autoregression parameters and thus remains insensitive to potential hidden linear trends or other structured disturbances, which frequently occur in economic, hydrological, and related applications.

Distributional shifts arise when the target deployment environment differs from the source environment that generated the training data. Robust learning frameworks such as Distributionally Robust Optimization (DRO) and Robust Satisficing (RS) aim to address this challenge, yet their finite-sample guarantees under such shifts, and their systematic comparison, remain underexplored: existing analyses typically establish guarantees either in the source environment or for adversarial worst-case performance over an ambiguity set. This paper instead studies generalization error in the target environment—the excess loss under the shifted target distribution. Our contributions are threefold. First, we derive finite-sample generalization error bounds in the shifted target environment for both DRO and RS. These bounds explicitly characterize the trade-off between reduced sensitivity to shift and the regularization penalty induced by each method's robustness hyperparameter, and they avoid the curse of dimensionality associated with Wasserstein empirical concentration. Second, when partial shift information such as shift magnitude or direction is available, we propose information-directed hyperparameter calibrations and compare the two methods given the same information. Under these calibrations, and in the partial-information regimes we study, DRO and RS exhibit complementary theoretical and empirical behavior. Finally, we apply the framework to a network lot-sizing problem, using it to interpret how robust policies respond to positive shifts in the demand distribution. Together, these results fill a gap in understanding the statistical properties of robust learning methods under distributional shifts and provide a principled basis for comparing DRO and RS.

We propose a robust barycenter for distribution-valued data by incorporating the Huber loss directly into the optimal transport cost. In contrast to metric-space Huber means, which apply the Huber loss to the Wasserstein distance after optimization, our construction acts on individual transport displacements, preserving quadratic behavior locally while limiting the influence of large displacements. The resulting Huber-Wasserstein barycenters form a natural interpolation between Wasserstein means and $L^1$-type Wasserstein medians. We establish the analytical and statistical foundations of this construction. For optimal transport with Huber loss, we prove regularity and uniqueness properties of dual potentials, existence of optimal transport maps, and stability as the Huber parameter varies. For the associated barycenter problem, we prove existence and characterization results, consistency of empirical plug-in estimators, and a finite-sample breakdown point essentially equal to $1/2$. In dimension one, we further derive the pointwise influence function and asymptotic distribution, quantify the associated robustness-efficiency trade-off, and show that displacement-wise Huberization can retain first-order information that is lost by distance-based Huberization under localized shape contamination. Numerical experiments on contaminated distribution-valued data demonstrate the robustness of the proposed barycenters and illustrate their interpolation between mean- and median-like behavior.

Mediation analysis is a powerful tool to decompose treatment effects into direct and indirect components, enabling explanations of total treatment effects in a formal statistical framework. However, applying such analyses to settings with time-to-event outcomes, time-dependent mediators and confounders remains challenging. Existing methods are statistically complex, computationally intensive, and rarely available in user-friendly software. The difference method offers a simple alternative, but its performance in this setting has not been systematically evaluated. We conducted a simulation study and real-world data analysis to fill this gap in the literature. Using Cox proportional hazards, Aalen additive hazards, and accelerated failure time (AFT) models with time-varying covariates, the difference method was compared across different data generation processes with time-dependent mediators and confounders, focusing on bias in estimated indirect effects. The parametric mediational g-formula was used as benchmark comparator. If correctly specified, the Aalen model based difference method produced unbiased estimates in the absence of a time-dependent confounder that was directly caused by the treatment. Similar results were obtained when using the Cox model based difference with rare outcomes, but not with common outcomes. The AFT model based difference method was biased in almost all scenarios, due to collapsibility issues. Only the parametric mediational g-formula was unbiased in all scenarios. In contrast to specialized methods, the difference method requires additional, often unrealistic, assumptions, such as the absence of a direct causal relationship of the treatment on time-dependent confounders. If those assumptions hold, however, it may be used as a simple and efficient alternative.

Machine learning theory studies learning procedures through mathematical setups in which the data model, training protocol, oracle access, loss, metric, and randomness define the phenomenon that a theorem is meant to explain. Solving an open problem therefore requires the problem formulation, theorem target, and proof mechanism to be developed in concert. Researchers formulate hypotheses, test them through preliminary theoretical or empirical analysis, and refine both assumptions and proofs. We investigate whether this process can be organized as an autonomous agentic workflow for ML theory research. We develop VALG, an agentic system that combines multi-level Verification, Adaptive formulation of Learning-theory problems, and Graph-structured proof development. Within each source-relative theorem branch, VALG maintains a fixed mathematical specification, checks the theorem-level composition of a typed proof-dependency graph, and constructs and reviews local proofs in dependency order. When a proof attempt fails, VALG identifies whether the obstruction lies in a derivation, the proof structure, or the theorem formulation and routes the next attempt accordingly. Formulation-level obstructions initiate an explicitly related variant or relaxation, preserving the mathematical relation between the resulting theorem and the source problem. We evaluate VALG on nine subproblems from five COLT 2026 open problems. Two runs produce internally finalized theorem candidates that match the scope of their source briefs; the remaining seven yield restricted-method results, special cases, or conditional theorems. These case studies show how VALG keeps source-scope matches, relaxations, conditional results, and blocked attempts mathematically distinct. VALG is open source at https://github.com/DechenZhang/VALG-ML-Theory-Agent.

This paper presents WIRED, an R package algorithm for joint probabilistic forecasting of multiple related time series. WIRED combines a library of simple marginal predictive distributions, CRPS-based adaptive mixture weights, and a Gaussian or Student t copula for cross-series simulation. We evaluate the implementation in a benchmark with four synthetic data-generating processes (DGPs), three forecast horizons, 30 replicates per DGP-horizon pair, nine ablations and external baselines, and a rolling-origin study on the built-in EuStockMarkets data. The central contribution is architectural and diagnostic. WIRED separates adaptive marginal expert aggregation from dependence reconstruction; the benchmark supports explicit dependence modeling, but shows that the current CRPS-extrapolated softmax weighting is not yet robust enough to dominate simpler bootstrap or equal-weight alternatives. The paper therefore identifies a working layer of the design, a bottleneck in the marginal aggregation layer, and a concrete research path for more regularized probabilistic ensemble construction.

In this paper, we study online statistical inference for quantile-based distributional policy evaluation, which aims to characterize the distribution of discounted cumulative rewards under a given policy. To obtain a finite-dimensional representation of the return distribution, we consider the fixed point $\bm\eta_m$ of the quantile-projected distributional Bellman equation, parameterized by the quantile vector $\bm\theta_m$. We develop two online inference procedures for $\bm\theta_m$, covering synchronous and asynchronous sampling from a generative model as well as Markovian data. The first method uses the trajectory of quantile temporal-difference learning (QTD). We establish functional central limit theorems for the partial sums of QTD iterates and construct a recursively computable random-scaling statistic. The resulting confidence intervals have asymptotically exact coverage at deterministic horizons, with length of order $T^{-1/2}$ in probability. The second method constructs anytime-valid confidence intervals directly from samples. Exploiting the monotonicity of the quantile Bellman map, we develop a computationally explicit coordinate-shaving algorithm with finite-sample simultaneous coverage. Under positive limiting state-visit frequencies, the resulting confidence intervals have lengths of order $\sqrt{\log\log T/T}$, which we establish as classwise optimal through a matching lower bound. Together, these results provide complementary guarantees for online inference in quantile distributional reinforcement learning, combining asymptotically exact fixed-horizon inference with finite-sample anytime validity.

Bayesian inference for A/B testing is a family of prior and stopping-rule configurations with fundamentally different statistical properties, but it is often discussed as a single method, and no systematic overview exists. This paper organizes common configurations into a three-tier hierarchy: 1) posterior coherence with no error control, 2) false positive rates bounded under continuous monitoring via Bayes factor stopping, and 3) false discovery rate control and calibrated shrinkage via empirical Bayes. Many commercial platforms operate at the lowest tier by default. We show that Bayes factor stopping is near-optimal for a broad class of cost functions, including most proposed in the A/B testing literature; because the same rule also controls the false positive rate, the choice between a decision-theoretic and a frequentist formulation is largely one of parameterization. Furthermore, the empirical Bayes prior is the only path to the third tier, but winner-selected corpora, pooled programs, and heterogeneous metrics can each prevent calibration regardless of corpus size. Simulations against group-sequential and always-valid frequentist baselines show that flat-prior posterior stopping exactly reproduces naive peeking, that a well-calibrated empirical Bayes prior achieves the lowest estimation error, and that expected-loss stopping minimizes regret only when shipping a null-effect variant is nearly free. Error rates, estimation accuracy, and regret are all different risks, and the appropriate method follows from the risks an experimentation program needs to control, not the other way around.

In this contribution we study max-stable random fields on the rooted tree under shifts to descendant subtrees. Branch-Brown–Resnick stationarity is characterised through homogeneous spectral classes, punctured tail measures, and local spectral tail fields. For lognormal representers, it is equivalent to invariance of the variogram under addition of a common prefix. We give Gaussian, max-autoregressive, regenerative cascade, and free-group cluster constructions, and show that summability on countably branching trees need not satisfy a zero–one law. We also derive the associated branch-invariant extreme-value and Archimax copulas.

Electrocardiography (ECG), photoplethysmography (PPG), and phonocardiography (PCG) provide complementary views of the same cardiac cycle, yet existing cardiac foundation models are trained for a single sensing modality, leaving the shared physiology across sensors unexploited. We introduce CardioState-JEPA, a cardiac foundation model to learn a single shared representation jointly across ECG, PPG, and PCG, built on a physiology-aware joint-embedding predictive architecture. The model maps heterogeneous waveforms into a common token space, processes them with a single shared Transformer encoder, and learns by predicting masked latent cardiac states, placing the pretraining target on shared physiology rather than sensor-specific waveform appearance. To handle the temporal offsets between electrical, mechanical, and hemodynamic events, cross-modal prediction uses a learned delay aligner that matches signals at the corresponding cardiac time. Because synchronized multi-sensor recordings are scarce, CardioState-JEPA first learns within-modality structure from abundant unimodal data and then uses paired data to align modalities in latent cardiac time. Evaluated as a frozen encoder across 25 downstream tasks spanning ECG, PPG, and PCG, our encoder improves average PPG classification by 8.2 AUROC points, PCG murmur detection by 18.8 AUROC points, and ECG classification by 15.5 AUROC points over the best self-supervised signal baseline and matches or exceeds cardiac models trained with privileged clinical text or supervised labels on several ECG benchmarks. These results establish that heterogeneous cardiac signals can mutually supervise a single foundation model of cardiac physiology.