New preprints, ranked for me

We study recursive-design wild bootstrap inference for dynamic panel data models with unobserved common factors estimated by Common Correlated Effects. In the large N,T setting, the bootstrap reproduces the biased limiting distribution in pure autoregressive models, but fails to capture all bias and factor-estimation variance components in models with additional regressors, particularly under weak exogeneity. We trace this failure to holding regressors fixed across bootstrap replications. We propose to combine bootstrap procedure with available bias-correction methods to conduct adjusted inference. Monte Carlo evidence shows substantial improvements over conventional strategies of using bias-correction paired with cross-sectional bootstrap methods.

In observational time series, statistical inference for dynamic causal effects of a one-time intervention across horizons is complicated by high-dimensional observed pre-treatment information, unmeasured confounding, and serial dependence. To address these challenges, we develop a semiparametric framework for inference from a single serially dependent time series, integrating debiased machine learning with instrumental variables through buffered block cross-fitting. Under geometric beta-mixing, we derive non-asymptotic bounds on estimation error, asymptotic normality at each fixed horizon, and feasible inference that accommodates serial dependence. We further show how learner-specific prediction guarantees under temporal dependence can be used to verify the nuisance-rate conditions required for orthogonal inference. In a monetary-policy application with 468 months and 1464 lagged FRED-MD controls, we show an instrumented policy tightening lowers housing starts at medium horizons, with sensitivity analyses that support the finding.

Information systems (IS) researchers increasingly use quasi-experimental methods such as difference-in-differences (DiD) and instrumental variables (IV) to recover causal effects from observational panel data. Power calculations that justify these designs assume i.i.d. errors, but the deeper problem is what even a cluster-robust calculator cannot see. We report a Monte Carlo study over 9837 parameter conditions (approx 9.8 million datasets) and decompose the planned-versus-achieved power gap. The serial-correlation component is recoverable by an AR(1)-aware calculator when rho is known, and partially when rho must be estimated from short pre-periods, but panel attrition, staggered-adoption bias, and parallel-trends pretesting are captured by no closed-form formula; exogenous attrition alone costs approx 8 to 11 percentage points at the few-hundred-to-thousand sample sizes IS studies use. Treatment-correlated, outcome-dependent attrition instead induces bias, not just power loss. For IV, holding first-stage F fixed, larger N neither raises power nor curbs exclusion bias, though with a fixed instrument more data does sharpen the first stage, so identification rests on instrument strength, not sample size.

We propose a nonparametric test for unobserved treatment effect heterogeneity in regression discontinuity designs. Under the null of no unobserved heterogeneity, a transformed outcome that imputes treated potential outcomes for untreated units must have a continuous conditional distribution at the cutoff. We convert this implication into an integrated conditional-moment restriction using characteristic functions, thereby allowing the conditional local average treatment effect to be an unrestricted function of covariates. We derive the asymptotic distribution of the test statistics via a $U$-process and establish the validity of a multiplier bootstrap procedure for calculating critical values. Monte Carlo experiments show well-controlled size and increasing power. Two empirical applications illustrate how the test distinguishes between heterogeneity explained by observables and that explained by unobserved factors.

This paper develops a general two-point dependent wild bootstrap (DWB) for weakly dependent estimating equations. Its key feature is that the two-point marginal distribution and the latent serial dependence specification can be chosen separately. The construction combines a normalized two-point distribution with a stationary latent Gaussian process via a Gaussian copula transformation, includes dependent Rademacher and Mammen multipliers, and nests the classical iid two-point wild bootstrap as the serially independent case. The induced multiplier autocovariances determine the lag weights in a corresponding heteroskedasticity- and autocorrelation-consistent (HAC) covariance estimator, which coincides exactly with the conditional covariance of the bootstrap estimating-equation sum. We establish first-order bootstrap validity for asymptotically linear estimators by showing that the matched-HAC estimator consistently estimates the long-run covariance and that the bootstrap estimating-equation sum converges conditionally to the same Gaussian limit as its original-sample counterpart, yielding valid HAC-studentized $z$-tests and the corresponding Wald and Lagrange multiplier tests. Monte Carlo experiments in nonlinear generalized method of moments (GMM) and linear regression show that Rademacher DWB generally provides more accurate finite-sample size control for $z$-tests than the Mammen and Gaussian DWB. A GMM application to a nonlinear short-rate mean-reversion model illustrates the practical relevance of the proposed two-point DWB.

A non-significant coefficient in a dynamic panel model need not imply the absence of a relationship. It may instead reflect heterogeneous effects averaged toward zero, reciprocal dynamics overlooked by a recursive specification, or relationships masked by the omission of correlated covariates. Standard linear estimators cannot distinguish among these possibilities. We develop an inferential workflow for political panel data that resolves this ambiguity by combining flexible autoregressive estimation, forecast-necessity testing, functional characterization, and same-data linear benchmarking. The workflow first identifies relationships required for out-of-sample prediction, then characterizes their functional form across political contexts, and finally, distinguishes differences arising from estimator flexibility from those due to model specification. Applied to the causal sequence model of democratization on the V-Dem panel of 113 countries, the workflow reproduces the model's central finding - the protective belt of civil society, the rule of law, and institutionalized parties - while recovering reciprocal relationships from democracy to its institutional supports that a linear model cannot detect. Most importantly, three weak published direct effects, of which two are null, and one is marginally significant, receive three different diagnoses: one dissolves under the full specification, one reflects heterogeneous effects averaged toward zero, and one was masked by the reduced variable set. The workflow corrects the published record in both directions, removing one relationship and recovering two. More broadly, the workflow provides a framework for evaluating dynamic political theories under a model class capable of representing nonlinear and reciprocal mechanisms while preserving relationship-level interpretation and explicit inferential standards.

It is standard practice to include covariates in regression discontinuity designs (RDDs) and regression kink designs (RKDs), but the theoretical justification for doing so does not generally extend beyond linear estimands. This paper proposes a novel entropy balancing reweighting approach for covariate adjustment within a general framework of RDDs and RKDs. While conventional regression-based covariate adjustment methods generally fail to deliver consistent estimation for nonlinear estimands such as quantile treatment effects, our reweighting approach achieves consistency while improving efficiency. Moreover, even in settings where the regression-based covariate adjustment method already improves efficiency, our approach can deliver additional efficiency gains. Simulation studies corroborate these theoretical findings. We present an empirical application in which our covariate adjustment yields statistically significant results that would not be obtained without covariate adjustment.

The tails of macroeconomic outcomes can respond differently from the centre of their distribution: shocks with modest effects on median growth or inflation can shift downside growth or upside inflation risk. We develop a threshold stochastic-volatility-in-mean VAR with regime-dependent leverage to study their structural drivers. The model allows endogenous interactions between outcomes and volatility, contemporaneous level-volatility dependence, and regime-specific propagation. In nearly 150 years of U.S. data, predictive model selection supports three inflation-defined regimes. We identify business-cycle, financial, macroeconomic-uncertainty, and financial-uncertainty shocks and decompose their contributions to growth- and inflation-at-risk. The structural composition of tail risk differs from that of the predictive median. Business-cycle shocks dominate the median response of GNP growth but account for a substantially smaller share of growth-at-risk. Macroeconomic uncertainty makes a material contribution to both growth- and inflation-at-risk, with its share of growth-at-risk increasing with the magnitude of a positive macroeconomic-uncertainty impulse, despite its limited role at the median. In high-inflation states, the contribution of financial uncertainty to inflation-at-risk rises with the magnitude of positive financial-uncertainty impulses.

Gaussian processes (GPs) provide a flexible nonparametric framework for modeling time series through appropriately chosen kernel functions. This chapter introduces the basic formulation of Gaussian processes, commonly used kernels, GP regression, hyperparameter estimation, and model evaluation using in-sample and out-of-sample criteria. Applications to stationary, quasi-periodic, and seasonal time series illustrate how individual and composite kernels can represent different forms of temporal variation. Additive kernels also provide interpretable decompositions into latent components such as trend, smooth local variation, and seasonality, while product kernels allow more complex dependence structures to be constructed. Finally, Gaussian process state-space models (GP-SSMs) are briefly introduced, and a nonlinear example demonstrates how a GP transition model can be combined with particle filtering and smoothing for latent-state estimation.

We develop a framework for identification, computation, and inference in econometric models with a linear-in-measures representation. These models express maintained restrictions as moment conditions linear in the joint probability measure of observed and latent inputs, and map that measure linearly to the distribution of outputs, even with nonlinear outcome equations. We construct an adversarial discrepancy function whose zeros characterize the identified set for structural and counterfactual parameters. With finite output support, finite linear programs compute the discrepancy function or provide certified bounds even when latent inputs have infinite support, and a penalized bootstrap yields confidence sets with uniform per-point coverage. We apply the framework to two open cases in binary choice panels with fixed effects and discrete covariates: sequential exogeneity with unspecified conditional marginal error distributions, and known conditional marginal error distributions with unrestricted serial dependence. In an entry game with multiple equilibria, the framework recovers the known sharp identification region.

Consider an optimal dynamic treatment regime, $g^{\textbf{opt}}$ correctly identified from a large, perfectly executed sequentially randomized experiment. Even when the experimental results are generalizable to a future target population, there is no guarantee that $g^{\textbf{opt}}$ outperforms human decision-makers; human experts can do better than $g^{\textbf{opt}}$ whenever they have access to relevant information beyond the covariates recorded in the experiment. Motivated by this observation, we derive results on a new class of regimes called initiation regimes, which generalize existing results on superoptimal regimes. These regimes follow human decision-makers up to the point where it becomes more beneficial to initiate a sequential optimal regime, and are guaranteed to outperform both purely human and purely algorithmic decision rules, e.g., based on reinforcement learning algorithms. Furthermore, we present modified experimental designs that identify the best initiation regimes, show how the best initiation regime can be identified from classical observational data under explicit assumptions, and give estimation and statistical inference methodology for these regimes. To illustrate the practical utility of the methods, we consider initiation regimes in a case study on treatment of lower back pain.

Many models, such as fixed-effect models for panel or network data, are hard to estimate because they feature nuisance parameters that are both numerous and estimated imprecisely. This, in general, causes an incidental-parameter problem in the estimator of the parameters of interest. The problem can be alleviated by working with an estimating equation whose expectation is insensitive to the value of the nuisance parameters. We discuss and contrast three notions of insensitivity, also called orthogonality, in the context of likelihood models: Neyman orthogonality, Neyman orthogonality to order q, and full orthogonality. Orthogonal moments are obtained by projecting the estimating equation on nested subspaces, which are spanned by, respectively, the scores of the nuisance parameters, the first q derivatives of the likelihood ratio with respect to the nuisance parameters, and all likelihood ratios of the model. We give explicit constructions in binary-choice, count-data, and nonlinear regression models.

We consider $K$-dimensional Bayesian vector autoregressions (BVARs) with Cholesky stochastic volatility (SV), in which the innovation covariance matrix is a lower-triangular linear transform of $K$ independent univariate SV processes. Such models are widely used in empirical macroeconomics to capture time-varying uncertainty and improve forecast accuracy, but the cost of posterior simulation is the binding constraint on the size of the system, and a major bottleneck for empirical work. We introduce a Markov chain Monte Carlo (MCMC) kernel that mixes better than existing samplers at the same computational complexity. It rests on a reparametrisation that makes the $K$ volatility trajectories conditionally independent, and on Particle Gibbs to update each trajectory. The kernel targets the exact posterior, rather than an approximation of it, and in an application with $K = 15$ it raises the mean effective sample size per second, relative to the benchmark corrected triangular algorithm, by a factor of approximately $14$ for the VAR coefficients and $3.4$ for the volatilities. We also introduce a a Sequential Monte Carlo squared (\smcsq{}) sampler, which used our MCMC kernel as a building block, and which delivers at every $t$ the one-step-ahead predictive density and the marginal likelihood of the data up to $t$, and hence on-line forecasting and model choice. To our knowledge, this is the first algorithm that delivers sequential marginal likelihoods for Cholesky-SV BVARs with static contemporaneous coefficients and a non-conjugate prior. We illustrate both on US monthly macroeconomic data, with $K=15$ for posterior inference and $K=6$ for the sequential sampler and model choice.

Many studies expose units to several binary treatments at once and record several correlated outcomes, yet analysts usually estimate one treatment's average effect on one outcome at a time -- discarding how treatments interact and how their effects vary across units. We introduce the factorial multivariate Bayesian causal forest, a Bayesian nonparametric model that decomposes the factorial response surface over the treatment lattice into a prognostic sum-of-trees plus one sum-of-trees per main or interaction effect, each with correlated multivariate leaf parameters, estimated jointly with coherent uncertainty. A single indicator-weighted kernel samples every component -- the prognostic term is the special case whose indicator is unity -- so the two-treatment model, the single-treatment multivariate causal forest, and a general order-r truncation are one and the same sampler. We give the ANOVA/M\"obius identification of each estimand, an efficient Rcpp engine, and an interpretability layer using value-suppressing uncertainty maps. In simulations the average-effect estimators are unbiased with nominal coverage, consistent, and robust under misspecification, failing only under unmeasured confounding, which we flag. We illustrate the method on clinical, agricultural and economic data and on a deeper application to the NHANES survey. Software is provided as an R package.

We propose a Bayesian framework for uncertainty quantification from the perspective that the working model is mis-specified in settings of a multilevel data-generating process. We focus on settings in which the mis-specification fails to match the functional form of the mean structure, and discuss Bayesian estimation of target parameters under dependence induced by a mismatch between working and data-generating models. The proposal represents a Bayesian semi-parametric procedure aimed at estimating population-level parameters while accounting for cluster- and unit-level variation in the estimating function. The proposal extends the regular Bayesian bootstrap to account for cluster- and unit-level variation using multilevel weights from an enriched Dirichlet model. Simulation studies indicate that the proposed approach has good frequentist properties when the data-generating process and the proposed model induce a partially exchangeable sequence associated with the unknown quantity of interest. Applications to radon (Gelman and Hill, 2007), Programme for International Student Assessment 2022 (OECD, 2023), and tuberculosis (Nobre et al., 2023) datasets are presented for illustrative purposes. The results demonstrate that the proposed method is competitive with variations of multilevel models, with major differences observed in the range of credible intervals, which are justified by the nonparametric assumptions underlying the proposed method.

I study treatment effect estimation when treatment events have persistent effects and can be experienced more than once. Natural disasters, job loss and health shocks are examples of such treatments. I show that the effect of a total treatment trajectory can be recovered under assumptions similar to those commonly invoked in single-event settings using suitably flexible TWFE models. Decomposing the total trajectory effect into portions attributable to distinct event occurrences, however, requires further assumptions. I propose an assumption similar to conditional parallel trends, imposing it on the growth of event-specific effects rather than on untreated outcomes. Combined with a linear-in-parameters model of effect growth, this assumption enables a sequential imputation estimator that consistently estimates the dynamic effects of each event occurrence and that can accommodate heterogeneity in effects according to observable event attributes, such as intensity. I demonstrate that several intuitive TWFE models fail to recover interpretable treatment effect parameters in the multi-event setting and illustrate the sequential imputation estimator's favourable performance using Monte Carlo simulations.

Policy reforms are sometimes accompanied by detailed individual-level data in the implementing jurisdiction, while only aggregate outcomes are available for potential comparison jurisdictions. This article develops an identification framework for heterogeneous policy effects when individual-level data are unavailable for the comparison units. The framework combines treatment-effect contrasts from difference-in-differences comparisons within the treated jurisdiction with a compatible population-average effect identified from aggregate panel data. Identification requires relative parallel trends within the treated jurisdiction together with the assumptions needed to identify the population-average effect from the aggregate panel. The within-jurisdiction component can be estimated from repeated cross-sections with a single pretreatment period.

Linear bandits model sequential decision-making problems with noisy rewards that are linear in the decision variable, where an agent must simultaneously learn about an unknown parameter that governs the mean rewards, while maximizing (expected) rewards over time. Two prominent algorithmic families--upper confidence bound (UCB) and Thompson sampling (TS)--achieve a balance of exploration (to estimate said parameter) and exploitation (utilization of knowledge about it) across time. The quality of estimation of that parameter depends on the eigenvalues of a design matrix. In this paper, we begin by showing that if the inference quality obtained from exploration, encoded in the minimum eigenvalue of the design matrix, grows $\gtrsim \sqrt{t}$ with time $t$, while actions remain sufficiently concentrated for exploitation, then an algorithm produces optimal high-probability $\mathcal{O}(\sqrt{T}\log T)$-regret rate over a time-horizon $T$ for spherical action sets. This analysis is algorithm-agnostic and follows an alternative route to the classical optimism-based elliptical-potential argument for regret analysis. Then, we illustrate that variants of UCB and TS satisfy the inference and concentration properties and in turn, enjoy optimal regret rate. In effect, our results provide a modular framework that can be used to analyze linear bandit algorithms and explicitly connect quality of parameter estimation to optimal regret accumulation.

Dynamic mode decomposition (DMD) is a data-driven method that computes the best linear approximation of the underlying dynamical system and decomposes the dynamics into a superposition of characteristic spatiotemporal patterns. Originally introduced by the fluid dynamics community, DMD and its extensions have found widespread use in many other research areas such as molecular dynamics, climate science, engineering, finance, and neuroscience. Applications include dimensionality reduction, forecasting, system identification, control, and spectral clustering. In order to apply DMD to partial differential equations, the spatial domain is typically first discretized using finite difference or finite element techniques, thus implicitly rendering the problem finite-dimensional. We extend projected and exact DMD to infinite-dimensional systems. Rather than estimating matrices from vector-valued observations, our DMD variants learn finite-rank operators from functional data such as observables, densities, or wavefunctions. We show that conventional DMD algorithms can be regarded as special cases of their functional DMD counterparts. All results will be illustrated with the aid of guiding examples. We focus in particular on Koopman, Perron-Frobenius, and Koopman-von Neumann operators associated with graphons, ordinary differential equations, and stochastic differential equations.

Price stability remains a pillar in monetary policy practices and carries a special importance within monetary unions. Mainstream economics tried to leverage price stability using price indices and several metrics to shed light on specific dynamics and optimal macroeconomic levels. The wide availability of data led researchers to consider the study of systems using Random Matrix Theory, based on inner correlation patterns. This aims to enhance the multivariate analysis by removing noisy patterns from the signal and improve data quality for further inferences. This work considers the collection of monthly inflation indices in the Eurozone as a \textit{system} of prices to analyze its eigenvalues' statistical and asymptotic properties and uncover inner country-level insights. Results confirm the system cannot assumed to be randomly generated, and the data exhibit noise-dominated patterns, due to small and persistent variations at the country-level. The latter make the inter-country correlations more dynamic and the separation of the signal from the noise quiet difficult. Findings identified two countries as distorting inflation dynamics besides three other distinct, regional-based groups of countries. Variability sources might stem from economic episodes fueling inflation spikes in some countries, as well as methodological aspects used to ensure data quality and representativeness in the European Union. Despite being complex, the system demonstrates a certain stability, in terms of self-organization; while large monthly fluctuations cannot be considered as rare events, but part of the data-generating process.

Automatic relevance determination (ARD) is the standard device for input selection in Gaussian process (GP) regression. By giving the covariance kernel a separate lengthscale for every input and learning those lengthscales by maximizing the marginal likelihood, ARD lets the data decide which coordinates matter: irrelevant inputs receive very large lengthscales and are effectively switched off. We trace this mechanism to the Bayesian Occam's razor embodied in the marginal likelihood, derive the gradient through which it prunes inputs, and emphasize that ARD delivers effective rather than exact sparsity. We review the algorithms used in practice and the rules that turn lengthscales into selections, and we survey the asymptotic theory, distinguishing the fixed-domain identifiability obstruction on the lengthscales from the high-dimensional selection-consistency guarantees recently established for hierarchical GP priors, and noting what remains open for plain ARD. We compare ARD with spike-and-slab priors, sparse axis-aligned and global-local shrinkage priors including the Bayesian lasso and horseshoe, penalized likelihood kriging, sensitivity and projection criteria, and additive kernels. We argue that ARD endures because of its seamless integration with kernel learning, universal software support, and low cost, and we close with its limitations and remedies.

Posterior sampling for reinforcement learning (PSRL) is one of the simplest and most effective exploration methods, but a basic question has remained open: does unmodified PSRL achieve minimax regret without structural assumptions on the prior? We answer yes. Exact vanilla PSRL is minimax optimal in leading-order Bayesian regret under arbitrary correlated priors. The difficulty is that a posterior-sampled transition model is coupled with its own continuation value. We overcome this with a common empirical transition reference that isolates the resulting value mismatch and a Bellman-based variance argument that controls it without an extra leading-order state-space factor. For finite-horizon, time-inhomogeneous tabular MDPs with unknown stochastic rewards, this yields the minimax $\widetilde{O}(\sqrt{SAH^3K})$ regret rate under arbitrary joint priors over rewards and transitions. The same proof principle gives the minimax $\widetilde{O}(d\sqrt{H^3K})$ rate for linear-mixture MDPs under arbitrary joint parameter priors.

We study fully online fixed-confidence causal discovery without any historical observational data. Starting from zero samples, the learner sequentially selects interventions to recover both the causal DAG and its edge weights under a linear-Gaussian structural equation model. We establish an instance-dependent lower bound for any $(\epsilon,\delta)$-correct algorithm and propose \textsc{FOCUS}, which adaptively allocates interventions through an online max--min game. A key contribution is a computable concentration inequality for the accumulated KL divergence involving causal parameters shared across interventions. We prove that \textsc{FOCUS} is $(\epsilon,\delta)$-correct and that its expected stopping time matches the lower bound in its $\Theta(\log(1/\delta))$ dependence up to an instance-dependent constant. Experiments demonstrate improved structure and edge-weight recovery and confirm the predicted stopping-time trend. Our codes are available on https://anonymous.4open.science/r/FOCUS_code-76E5

Consider testing a finite composite null $\cP={P_1,\ldots,P_L}$ against a simple alternative $Q$ using $n$ i.i.d. observations. Let $\ell_n$ denote the supremum of the expected log e-value over e-variables that are exact under every $P_i$ and pivotal across the nulls. Zhang, Ramdas, and Wang (2024) showed that $\ell_n$ is superadditive and asked whether $\ell_n/n$ converges to the upper bound $\min_i D(Q|P_i)$ under their joint atomlessness and absolute continuity assumptions. We prove that it does. The proof draws on techniques developed by Zhang, Ramdas, and Wang (2024) and by Farooq, Fritz, Haapasalo, and Tomamichel (2024).

In this work, we study the one-dimensional regression problem under random design and Gaussian errors. Our framework is very general: we make no prior assumptions about the design (which may be nonstationary and exhibit short or long-range dependence), nor do we assume that the errors are homoscedastic. We examine in detail the cases where the error process exhibits short or long-range dependence. We adopt a least-squares penalized strategy using piecewise polynomials to estimate the regression function, following the framework of Baron, Birg{\'e} and Massart [1999]. We derive explicit penalties, up to calibration constants, to obtain adaptive estimators for which we establish risk bounds. Since these penalties depend on the dependence properties of the error process, which are unknown in practice, we propose several adaptations of the dimension jump calibration algorithm to make our procedures fully data-driven.

This paper presents an end-to-end causal machine learning (ML) pipeline designed for real-world applications with continuous treatments. The proposed framework consists of six sequential steps: dimensionality reduction, causal identification, positivity assumption violation handling, estimation, refutation and evaluation, and policy optimization. We introduce practical contributions not currently available in existing causal ML toolkits, specifically: (1) a method for detecting and quantifying positivity violations in continuous treatment settings (2) a novel, scalable two-stage dimensionality reduction framework tailored for causal inference with high-dimensional data; (3) the adaptation of sensitivity analysis and estimation methods originally designed for binary treatments to the continuous treatment space and (4) an end-to-end integration of these components into a modular, reproducible workflow. These innovations address real-world challenges in causal inference that are often not covered in theoretical frameworks but frequently encountered in industrial applications. The methodology is validated with a synthetic dataset inspired in a real-world financial debt collection use case, however its design can be applied to analogous problems across different industries. Results demonstrate that the proposed methodology offers a more computationally efficient approach and produces less biased estimates compared to standard methods for problems with continuous treatment and high-dimensional data. A fully functional GitHub repository with documented code and numbered notebooks is made available ensuring reproducibility and practical implementation. The pipeline presented is intended to contribute to closing the gap between academic approaches and practical application in industry contexts where causal ML can be highly beneficial such as the financial sector.

Martingale posteriors and related predictive resampling methods replace the likelihood--prior pair used within Bayesian inference with a predictive model for future observations. These methods are simple to implement and increasingly popular due to their computational efficiency, but little is known about their ability to accurately quantify uncertainty. In this work, we study the concentration and calibration properties of the martingale posterior for general functionals, and show that credible sets can systematically undercover if the predictive algorithms are not carefully tuned. We propose a simple remedy: the bagged martingale posterior. Rather than starting every predictive path from the observed sample, we start the paths from random bootstrap resamples of the data and then amalgamate the resulting draws. Critically, this approach incurs no additional {simulation} cost compared to standard predictive resampling, and delivers conservatively calibrated credible sets. This scalability enables application to a range of challenging examples, including sparse high-dimensional regression and nonparametric conditional quantile models.

We consider the construction of confidence intervals for population means when observations are subject to missingness. To accommodate more general missingness mechanisms than missing completely at random (MCAR), we adopt a reparametrised version of the realisable contamination model of Ma et al. (2026), which is a mixture of an MCAR version and a missing not at random version of the same base distribution $P$. We characterise the minimax length of confidence intervals that can adapt to potentially unknown parameters of the model, including the contamination fraction, for Gaussian base distributions and for nonparametric classes satisfying certain tail or symmetry assumptions. In all of these settings, we provide explicit constructions of simple, practical and finite-sample valid adaptive confidence intervals that attain the corresponding minimax rates. Finally, we provide implications of our results for causal inference.

In this paper, we formulate a general differencing framework for variance estimation across a range of settings. We demonstrate that conventional difference-based noise variance estimators cannot achieve the desired bias-correcting power in nonparametric regression with repeated measurements. A new high-order bias-corrected differencing scheme, adapted to repeated measurements, is proposed by interlacing inter-group and intra-group differencing. The theoretical properties of the new sequences and estimators are studied. Our proposals are particularly efficient in finite samples and under high signal-to-noise ratio scenarios, where asymptotic convergence has not yet fully taken effect, due to their strong bias-correcting power.

Bayesian Federated Inference (BFI) estimates statistical models from multicenter data when individual-level observations cannot be combined across centers. We present \pkg{BFI}, an \proglang{R} package that implements this methodology for Gaussian, binomial logistic, and survival regression models. Each center performs a Bayesian maximum a posteriori analysis and sends only parameter estimates and curvature information to a central server, where these summaries are combined to approximate the analysis of the combined data. The package supports prior specification, structured forms of between-center heterogeneity, several parametric and flexible baseline-hazard models for survival analysis, and treatment-effect estimation for observational and randomized studies. We describe the software design and the information exchanged between centers, give reproducible workflows for these analyses, and compare the federated results with pooled-data analyses for Gaussian, logistic, and survival models.

Panel conditioning, the causal effect of prior survey participation on responses, can vary with tenure. Under an additive model of cell means in period, entry cohort, and tenure, we characterize which features of the conditioning path a staggered panel identifies on its observed support, and how the unidentified component affects common panel estimators. The identified set of the path is an affine translate of the tenure projection of the cell design's kernel, and a linear functional of the path is identified exactly when it annihilates that projection. It always contains an affine direction and, when the entry cohorts share a stride, periodic directions, which exhaust it under a connectivity condition on observed increments; second differences at that stride are then identified, and ordinary ones generally are not when the stride exceeds one. Under a recruitment condition, an interrupted schedule such as the four-eight-four rotation of the Current Population Survey (CPS) distinguishes a constant increment per interview from one per calendar month, which no equally spaced schedule can. We give support conditions for recovery under a plateau, entry-wave negative controls, or bounded cohort drift. A second set of results links identification to regression: two-way fixed effects absorb every unidentified direction, so the remaining conditioning bias is normalization-invariant and itself identified, and a two-way regression with tenure indicators corrects it under a residual-rank condition. When event time is aligned with tenure, conditioning shifts event-study coefficients by a known linear functional of the path, producing pre-trends without anticipation; bounds on identified curvature bound those shifts. Simulations verify the identities, a 19-wave Japanese panel illustrates the support calculations, and published CPS month-in-sample indices give a descriptive, not identifying, example.

We develop a general framework for model specification testing in stationary functional time series. The approach is based on an autoregressive approximation that represents a broad class of stationary functional processes through coefficient kernels whose dimension and autoregressive order may increase with the sample size. Different model assumptions induce different structural restrictions on these kernels, and our tests are constructed by measuring deviations from the corresponding restrictions. We illustrate this principle for three problems: testing a prescribed order of a functional autoregressive model, testing a functional autoregressive moving-average specification, and testing separability of autoregressive coefficient kernels. The resulting statistics are based on weighted $\mathcal{L}^2$-distances, and the critical values are obtained by a multiplier bootstrap. We establish a quantitative bootstrap approximation that is uniform over a class of weight functions and prove asymptotic validity and consistency of the proposed tests. The methodology allows for data-adaptive weighting and is illustrated by simulations and a data example.

Instead of having a single "yes" or "no" result from a test of the global null hypothesis that a function is increasing, we propose a multiple testing procedure of the function's increasingness at several points. If the global null is rejected, then multiple testing provides more information about why. If the global null is not rejected, then multiple testing can provide stronger evidence in favor of increasingness, by rejecting null hypotheses that the function is decreasing. Our approach uses high-level assumptions that apply to a broad class of causal and descriptive statistical models. By inverting the proposed multiple testing procedure that controls the familywise error rate, we also generate "inner" and "outer" confidence sets for the set of points at which the function is increasing. With high asymptotic probability, the inner confidence set is contained within the true set, whereas the outer confidence set contains the true set. We also improve power with stepdown and two-stage procedures. Simulation and empirical examples illustrate the new methodology, and all code is provided.

For decades, statistical inference has largely been developed one problem at a time. Given a scientific target, such as a treatment effect or a regression function, statisticians design a problem-specific estimator together with a procedure for quantifying its uncertainty. This paper proposes a different paradigm. We focus on a classical problem in statistical inference, confidence interval construction, and develop TabCon, an amortized inference system built on a tabular foundation model that produces confidence intervals for new datasets through a simple forward pass. The key methodological ingredients of TabCon are a sparse mixture-of-experts architecture and reinforcement-learning-based post-training that calibrate the resulting confidence intervals to a desired coverage level. Across a wide range of benchmark datasets, TabCon attains near-nominal coverage while producing short confidence intervals. At inference time, it also offers considerably greater computational efficiency, running 50 times faster than the classical bootstrap procedure, even when the latter uses only 50 bootstrap samples.

Time series forecasting requires extrapolating the dynamics of an observed process beyond the last available measurement. Yet recurrent forecasting models typically treat observations as inputs that directly control their latent dynamics. It leads to a regime change when these observations become unavailable at prediction time. Following a state-estimation perspective, we introduce Observer State-Space Models (OSSMs), a class of recurrent models that interprets the observed input time series as measurements of an underlying autonomous dynamical system. OSSMs explicitly separate latent-state propagation from measurement assimilation: a single transition governs the dynamics across both context and forecasting intervals, while available observations correct the estimated state through an observer. This formulation naturally exposes classical control-theoretic properties, including observability and convergence of the state estimation error. We further show that conventional and recent SSMs can be recovered as particular instances of our OSSM framework, thereby providing a unified interpretation of their recurrent dynamics and revealing modeling inconsistencies. We perform experiments across several benchmarks showing that OSSM achieves substantial improvements while maintaining the same parameter count and training setup as the corresponding SSM baseline. These results support a simple principle for recurrent forecasting: observations should correct the estimated latent state, rather than control the dynamics used to propagate it.

Prediction-powered inference (PPI) typically relies on point predictions on unlabeled data. When predictive distributions are available as in a wide range of applications, including predictions from LLMs, we introduce distribution-informed prediction-powered inference (DiPPI), a general framework for further improving statistical efficiency by using predictive distributions as auxiliary information in the spirit of PPI. We characterize the optimal use of this information through score calibration, derive the oracle efficiency for a finite-dimensional representation of the predictive distribution, and provide theoretical guarantees for positive learning with the cross-fitted DiPPI estimator. Through simulations and three real-data applications, we show that DiPPI achieves better efficiency than PPI methods based on point predictions. These gains arise when the predictive distribution contains score-relevant information that is partially lost in point predictions.

We derive the limiting distributions of the $M$-test family of unit root statistics in the nearly integrated nearly white noise (NINW) framework introduced by Nabeya and Perron (1994) in the case of an unknown linear time trend. In the case of known long run variance (LRV), the limiting distributions of the $M^{GLS}$ tests are contaminated by additional noise terms as a result of quasi differencing whereas these terms are less present in the $M^{OLS}$ limiting distributions, both of which display conservative properties under conventional critical values. Furthermore, we prove the Gaussian power envelope in the NINW model is asymptotically equivalent to the standard envelope of Elliott, Rothenberg, and Stock (1996), and that the oracle $M$-tests have inefficient power relative to this benchmark. We then derive the limiting distributions of the feasible statistics and show that the autoregressive estimate of the LRV commonly used overestimates the LRV, creating altered limiting distributions. Finally, finite sample simulations illustrate that, of the procedures considered, no uniformly satisfactory solution exists for handling a series with a large negative moving average coefficient.

The Baringhaus-Henze-Epps-Pulley (BHEP) tests for multivariate normality are affine-invariant goodness-of-fit tests based on a Gaussian-weighted $L^2$ distance between empirical and Gaussian characteristic functions. In 1990, Henze and Zirkler expressed the limiting null distribution through the eigenvalues of an integral operator on the standard Gaussian space. In 1997, Henze and Wagner obtained a simpler covariance kernel and raised the problem of calculating the eigenvalues of the resulting operator on a Gaussian-weighted space. Although subsequent work treated the univariate case and numerical approximations in a few low dimensions, the complete all-dimensional spectral problem remained open. This paper determines both complete spectra for every dimension $d \in \mathbb{N}$ and every smoothing parameter $\beta > 0$. The two operators are shown to have the forms $\mathcal{X}_{\beta,d}^*\mathcal{X}_{\beta,d}$ and $\mathcal{X}_{\beta,d}\mathcal{X}_{\beta,d}^*$ for the same Hilbert-Schmidt operator $\mathcal{X}_{\beta,d}$. Consequently, their nonzero eigenvalues agree, including multiplicities, while the null space of the Henze-Zirkler operator is identified exactly. The Gaussian integral operator in the Henze-Wagner decomposition is diagonalized by Mehler's formula, and rotational symmetry confines the finite-rank correction to the sectors associated with spherical harmonics of degrees $0$, $1$, and $2$. The degree-$1$ and degree-$2$ eigenvalues are characterized by scalar transcendental equations, and the radial eigenvalues by an explicit pole-safe Fredholm determinant. The paper establishes nonnegativity, multiplicities, eigenfunction reconstruction, completeness, the trace identity, and a complete characterization of all exceptional pole cases.

Automatic relevance determination (ARD), the default tool for variable selection in Gaussian-process (GP) regression, ranks inputs by inverse lengthscales -- which measure how fast a function varies, not how much an input contributes to prediction -- and offers no calibrated rule for deciding which inputs to keep. The prediction-centred alternative, the derivative sensitivity $\nu_j = \mathbb{E}[(\partial f/\partial x_j)^2]$, is available in closed form from a fitted GP, but turning it into a selection rule is harder than it looks: at a null input the estimator is a degenerate quadratic form, so Wald and Bernstein-von Mises cutoffs are anti-conservative, and the natural residual bootstrap is mis-scaled. We show that a studentized multiplier bootstrap of the GP derivative process repairs both, prove its validity through an invariance principle for quadratic forms, and obtain asymptotic family-wise and false-discovery-rate control across inputs. Over 100 replications the rule controls FDR wherever inputs are truly null, while uncalibrated derivative rankings breach the target by up to 2x and a Bernstein-von Mises cutoff by 2.2x; at matched FDR it loses no power; it holds under a Mat\'ern kernel and input correlation up to 0.99; on real data with planted and authentic null inputs it admits 5-12x fewer spurious inputs; it costs 5-18% of the GP fit; and a block-averaged variant retains validity at cost linear in $n$.

Rubin's variance ratio (VR) complements the standardized mean difference by detecting covariate imbalance in spread, but no purpose-built omnibus extension exists for more than two groups. We introduce FVR, a size-weighted quadratic combination of pairwise log-variance-ratios that generalizes VR to k groups, together with geometric-mean and maximum pairwise variance-ratio comparators. FVR has an exact relationship to Rubin's VR at k = 2 and a mathematical structure directly parallel to Cohen's f. In Monte Carlo simulations spanning four variance-driven bias mechanisms, k = 3, 4, 6, and sample sizes from 200 to 10,000, FVR and the geometric-mean statistic generally tracked downstream estimation bias at least as well as the maximum statistic. FVR retained substantially more signal when many groups shared a very small sample, while absolute-bias thresholds were unstable at N = 200. Across mechanisms, FVR < 0.10 was reasonably reassuring, values above approximately 0.30 generally indicated concern, and intermediate values were context-dependent. The statistic is implemented for arbitrary analysis weights in the Stata command varatio and illustrated using a multivalued-treatment application.

Competing-risks data do not identify latent marginal duration distributions or their dependence without additional restrictions. This paper asks whether restrictions introduced to restore identification themselves restrict the observable law. For a two-risk Archimedean model with categorical exclusion restrictions, we derive a necessary-and-sufficient observable characterization. A discrete single-crossing argument identifies the scalar copula parameter from cell-specific overall survival probabilities, while cause indicators recover the remaining allocation and generate additional specification restrictions. We construct an identification-robust, self-normalized quadratic statistic, invert it to obtain confidence sets, and use empty inverted sets as a conservative specification test. Simulations show that the cause indicator can be decisive: under the weakest contrast considered it converts a frequently uninformative survival-based confidence set into an informative joint set without loss of coverage, whereas excessive categorical contrast can eliminate causes from individual cells and make cause-specific recovery inadmissible. The results turn latent exclusion restrictions into refutable restrictions without estimating covariate derivatives.

Unobserved common causes are pervasive in real-world time series and can induce spurious associations that causal discovery methods mistake for direct edges. We propose LUCID (Learning Under Confounding for Inference and Discovery, a regime-adaptive deconfounding layer that first estimates the confounding regime from data using a Mar\v{c}enko--Pastur spectral router, then applies a deconfounding strategy matched to that regime. When the spectrum indicates pervasive factor confounding, LUCID attenuates factor-dominated variation and recovers contemporaneous (lag-$0$) structure from the resulting innovations, with edge selection calibrated against a data-driven edge-free null. Rather than being tied to a particular discovery algorithm, it can wrap existing discovery engines; we demonstrate consistent improvements across three such methods. On a diverse synthetic out-of-distribution benchmark spanning changes in confounder strength and sparsity, loading density, lag structure, volatility dynamics, edge heterogeneity, persistence, intermittency, and tail behavior, LUCID achieves the best family-weighted directed, lag-resolved graph $F_1$ ($0.60$), improving over the strongest baseline by $0.19$ absolute ($\approx\!46\%$ relative). Its advantage widens relative to looser lag-collapsed scoring, and remains robust under intermittent and heavy-tailed confounding. Code reproducing the method, the benchmark generators, and every reported experiment is available at https://github.com/bloomberg/causal-ts.

Staggered distributional difference-in-differences produces cohort-specific potential-outcome distributions, but applied work typically wants one overall quantile treatment effect. Two natural summaries---averaging cohort quantile treatment effects (QTTs) and mixing cohort distributions before inversion---use the same policy weights yet answer different target-population questions and can disagree even in sign. We derive the exact sharp interval for their gap conditional on cohort quantiles and weights, a globally sharp range-only envelope, and a locally sharp density-tilt representation. We develop joint smooth and mass-point-safe inference and show that neither estimator is uniformly more precise, even when the estimands coincide. In a same-object reconstruction of a public staggered-QTT application, holding data, identification, distributions, and weights fixed while changing only aggregation order reverses reported signs at several quantiles. Aggregation order is therefore part of the estimand and must be chosen before inversion.

Tool-using agents are continually updated with new interaction data. After each policy update, however, previously estimated action credits may become stale. Recomputing them from scratch can require many additional tool calls and environment interactions, making repeated updates increasingly expensive. We ask a simple question: when does historical action credit actually need to be updated? Our key observation is that a change in action value does not necessarily imply a change in the decision. Historical credit can still be useful as long as policy-induced drift is too small to overturn the existing action ranking. Building on this idea, we introduce pairwise branch sensitivity to capture how strongly a policy update affects the downstream regions that distinguish two candidate actions. We then derive a first-order anchored credit-transport estimator that updates historical credit using old interventional trajectories, and propose a Decision-Sufficient Credit Gate (DSC-Gate) that chooses whether to reuse, transport, or resample credit. Experiments show that branch sensitivity explains credit drift substantially better than global policy distance. With sufficient historical data, credit transport reduces estimation error, while its benefit to decision making is concentrated on updates that affect action-distinguishing branches. On a fully independent test set, DSC-Gate changes mean regret by only +0.00004 relative to a gap-based gate while reducing mean new tool steps from 472 to 286, a 39.4% reduction. We observe the same pattern after a real tool-agent parameter update. Overall, our results show that agents do not need to recompute action credit after every policy update: much of the historical evidence can be reused or cheaply corrected, reducing the additional interaction required to keep action decisions up to date.

Prediction with expert advice is a fundamental problem in online learning. When the time horizon $T$ is known in advance, the minimax cumulative regret over $n$ experts is asymptotically $\sqrt{\frac{T \ln n}{2}}$. This is achieved by the Multiplicative Weights Update algorithm with a learning rate tuned to $T$, and is known to be tight. If instead the regret bound is required to hold simultaneously at every time $t$, the best known guarantee has been $\sqrt{t \ln n}$---a factor of $\sqrt{2}$ worse---and it has remained unknown whether this factor of $\sqrt{2}$ is necessary. We show that it is not. We give an algorithm, requiring no knowledge of the horizon, whose cumulative regret satisfies $R_t \le \bigl(1 + O(\sqrt{\ln \ln n / \ln n})\bigr)\sqrt{t \ln n / 2}$ simultaneously for every $t \ge 1$.

This paper conducts an empirical mediation analysis to examine how U.S. stock-market sectors responded to the 2025 tariff-policy episode. Using S&P 500 stock data along with hundreds of firm-level financial variables, we study whether sector membership was associated with differences in stock returns and whether these differences were statistically related to firm financial characteristics. For each of five economically distinct tariff-policy windows, we use high-dimensional linear causal mediation models to estimate the direct sector effects. We use a linear regression model to estimate the total sector effects, and then estimate the indrect sector effects by using the difference between the total sector effects and the direct sector effects. To accommodate a large number of candidate mediators, we use a partially penalized least-squares procedure that regularizes financial-variable coefficients while unpenalizes coefficients of all sector indicators. Wald and F-type tests are used to examine whether indirect and direct effects are significant or not, respectively. The results show that tariff news did not produce one uniform market response. Sector effects vary across five economically distinct tariff-policy windows. The clearest findings are in the direct sector effects. Technology was the sector most negatively affected during the initial decline and escalation collapse, but it also rebounded most strongly during the April 9 relief rally. The results indicate that in addition to Technology, Energy and other trade-exposed sectors were particularly affected during the initial decline and escalation phases. The selected mediators provide evidence that growth, investment, liquidity, profitability, and cash-flow characteristics helped characterize firms associated with sector-level return differences.

Understanding why predictors can generalize despite interpolating noisy training data is a central puzzle in machine learning. Most work on such "benign overfitting" studies minimum-2-norm linear regression, reflecting the inductive bias of gradient descent. However, modern optimizers such as Adam and Muon use non-Euclidean update geometries, favoring solutions associated with other norms. Analyzing regression for non-Euclidean norms is substantially more difficult, with known results essentially limited to Gaussians. In this paper, we develop a method to analyze benign overfitting in linear regression for general norms and general (sub-Gaussian) distributions. As a special case, we prove that minimum-p-norm interpolation with p>1 can benignly overfit even for non-Gaussian distributions, under suitable conditions. Perhaps surprisingly, for the 1-norm, benign overfitting does not hold in general for well-behaved (but non-Gaussian) distributions, showing that existing positive 1-norm results rely crucially on Gaussianity. Our proof analyzes the geometry of the dual optimization problem, using concentration and central limit tools to show it is approximately Euclidean in many high-dimensional cases.

A central goal when designing treatment policies is often to "do no harm", that is, to avoid interventions that improve average outcomes while worsening outcomes for some individuals. A widely used notion for harm is the fraction of negatively affected (FNA), defined as the probability that an intervention decreases an individual's outcome. However, in many applications, treatments operate through mediators, and a single "total" FNA can obscure whether harm arises primarily through direct pathways or indirect (mediator-induced) pathways. In this work, we introduce a path-specific analogue of the FNA. For this, we disentangle total harm into direct and indirect harm in causal mediation settings. However, these quantities depend on joint distributions of potential outcomes that are not point-identified even in randomised controlled trials. As a remedy, we develop a novel partial identification framework for direct and indirect FNA. In our framework, we (i) derive sharp Makarov bounds for the FNA, and (ii) propose a semiparametrically efficient estimator with valid confidence intervals for these bounds under mild margin conditions. We demonstrate our framework across various numerical experiments. To the best of our knowledge, we are the first to study path-specific decomposition of causal harm and to develop an orthogonal inference framework for its analysis.

The proportion mediated (PM), defined as the indirect effect divided by the total effect, is intuitive when direct and indirect effects operate in the same direction. Under inconsistent mediation, however, opposing pathways partially cancel each other, causing PM to fall outside the unit interval or become unstable as the total effect approaches zero. We introduce two complementary effect-scale measures: the Mediated Magnitude Share (MMS), which quantifies the absolute magnitude of an indirect pathway relative to the combined absolute magnitudes of the direct and indirect pathways, and the degree of suppression, which quantifies the fraction of this combined magnitude that is cancelled by opposing pathways. We establish the bounds and limiting behavior of both measures in the two-pathway case, show that MMS reduces to PM under consistent mediation, and extend both measures to settings with multiple indirect pathways. For binary outcomes, we extend the framework to odds-ratio and risk-ratio scales by using logarithms to convert multiplicative effects into additive ones, enabling direct application of MMS and the degree of suppression. Together, these measures provide bounded and interpretable characterizations of mediation systems in which pathway directions are not uniform.

We investigate the least squares linear regression problem with random partial Discrete Fourier Transform (DFT) matrices, providing a rigorous analysis of the model's generalization error. By leveraging tools from random matrix theory, we derive exact non-asymptotic bounds for the risk of the Moore-Penrose estimator, which hold for finite-dimensional problems and reveal the precise dependence on key parameters such as the sample size, dimension, and noise variance. Then we obtain a characterization of the double descent phenomenon in the linear regression context, demonstrating how the risk evolves when the number of parameters $p$ and the number of samples $n$ tend to infinity, with $p/n$ fixed. The analysis relies on applications of the Stieltjes transform for random Fourier matrices, enabling a precise description of the spectral properties of these matrices and their impact on regression performance. To validate our theoretical findings, we present several numerical examples that illustrate the double descent curves. These simulations align closely with our derived bounds, confirming their predictive power in both under-parameterized and over-parameterized regimes.

An analyst should not assume that a design score improves a meta-analysis simply because it is labeled as robustness. This paper treats declining heterogeneity along a prespecified, outcome-separated design score as a falsifiable directional hypothesis. In an 18-study illustration, the constrained model lands at $\gamma = 0$ and becomes numerically identical to conventional random effects, while the unrestricted scale fit points in the opposite direction. Simulations that deliberately violate the restriction show the same graceful fallback under wrong-direction and U-shaped heterogeneity; a stepwise decreasing variance function retains an efficiency gain despite functional-form misspecification. DR-Meta formalizes the hypothesis as $\tau^2(\mathrm{DR}) = \tau_0^2 \exp(-\gamma \mathrm{DR})$, $\gamma \geq 0$, within an established location-scale framework. Exact nesting, conditional weight monotonicity, a fixed-effect upper bound, pseudo-true targets under mean misspecification, and score-rescaling behavior are derived with full proofs. Efficiency gains are negligible around an empirically anchored $\gamma$ near 1 and become material only under stronger gradients; the largest RMSE reduction, about 9.6 percent, occurs in a deliberately extreme condition. Plug-in intervals under-cover relative to random effects, and modified Knapp-Hartung adjustment improves coverage only modestly. Widening the $\gamma$ optimization bound changes $\gamma$ estimates much more than pooled-estimate RMSE. The practical contribution is therefore diagnostic rather than automatic efficiency: a nested directional restriction can be supported, contradicted, or collapse transparently to random effects, provided constrained, unrestricted, location, and sensitivity analyses are reported together.

Two-sample Mendelian randomization (MR) is a widely applied methodology in epidemiology. In two-sample MR, summary data (typically, regression coefficients and standard errors) quantifying the association between multiple genetic variants and the exposure and the outcome are used in an instrumental variable framework aimed at estimating the causal effect of the exposure on the outcome. Most two-sample MR methods were developed under data-generating models where the association of for each candidate genetic instrument with the exposure, as well as the causal effect of the exposure on the outcome, are constant in the additive scale. These assumptions are useful because they imply that, had all genetic variants been valid IVs, they would all estimate the same causal parameter - namely, the constant causal effect. We refer to this condition as summary-level homogeneity. However, these are rather strong homogeneity conditions which may raise concerns about the plausibility of these methods in practice. In this paper, we show that summary-level homogeneity is implied by the following conditions: the causal effect is additive linear, but not necessarily constant across, all strata of the population; and uncorrelatedness between heterogeneity in the causal effect and in the association between each genetic variant and the exposure. Under these conditions, typical two-sample MR methods can be interpreted as estimators of the average causal effect. These results clarify that point-identifying assumptions required for two-sample MR methods are weaker than previously anticipated, which contributes to their plausibility and interpretation in at least some practical applications.

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

A special semimartingale admits a unique decomposition $X=X_0+M+A$ into a local martingale $M$ and a predictable finite-variation part $A$. We consider the identification of $A$ when $X$ is observed only through repeated cross-sections. The estimand is then the projection of the sampled predictable compensator onto the observable feature filtration, namely the current state together with whatever randomness is shared across the population, so that at a fixed diffusion coefficient the marginal flow identifies the drift only up to a Markovian projection. If the drift is an affine functional of an observed lag window, the joint problem is a convex quadratic programme whose solution is the pseudo-panel regression of econometrics. Our principal concern is the case, which we believe not to have been treated before, in which the drift is the output of a hidden linear dynamical system whose dynamics are themselves to be identified from the marginals. The joint problem is then a bilinear quadratically constrained programme, which we solve to certified global optimality by spatial branch and bound; with unpenalised state disturbances and a drift basis growing with the grid it is NP-hard already in latent dimension one, by reduction from $\ell^1$ rank-one matrix approximation, whereas the complexity of the deterministic system at fixed latent dimension remains open. A block-coordinate decomposition offers a cheaper alternative. For the estimator itself, we obtain rates at a fixed mesh, separated into Monte-Carlo, estimation and grid contributions.

The SLGP package provides nonparametric estimation of fields of probability densities indexed by covariates, from unevenly sized samples with few or no replicates. A Spatial Logistic Gaussian Process applies a logistic density transformation to a finite-rank Gaussian process over a joint index-response domain, yielding a prior over conditional densities that adapts in location, shape, and modality. The contribution of this paper is the implementation: a finite-rank Random-Fourier-Feature representation that keeps the transformed process numerically tame, a log-concave likelihood that makes the fit a convex problem, and a family of quadrature shortcuts for the normalising integrals that turn a per-observation cost into a per-grid cost. Three estimation modes are available: MAP, Laplace, and full MCMC. We demonstrate the package throughout on the classical Fiji-Tonga earthquake catalogue, estimating how the distribution of hypocentre depths deforms across the region, and on a discrete-response variant.

This is mostly a non-technical note, reflecting the author's philosophy and personal views on time series analysis, with references limited to only a few representatives for brevity. The Threshold Principle, formally announced in Tong (1990), modernised time series analysis by introducing a collection of sub-systems to model complex nonlinear dynamics. We explore the conceptual architecture of the Threshold Principle drawing parallels and contrasts with artificial neural network in machine learning. We trace the influential adoption of Threshold Autoregression in econometrics. We examine a smooth extension, namely the smooth threshold autoregression introduced by Chan and Tong (1986) that was later popularized in the Econometric literature. We sound cautions to help econometric users to avoid misuse of this model. Furthermore, we examine how we can systematically apply the Threshold Principle to conditional variance to enable meaningful volatility classification. Finally, we mention some of the modern applications of the Threshold Principle to non-real-valued domains, underscoring its enduring half-century methodological significance as embodied in the threshold autoregression.

Off-policy evaluation is a foundational component of offline reinforcement learning, aiming to assess and optimize policy performance using pre-collected datasets. However, such datasets often suffer from pronounced challenges, including distribution shift, $Q$-value overestimation, and low sample utilization efficiency. To address these issues, this paper introduces a weighted Bellman residual minimization framework that incorporates density ratio weighting by effectively integrating expert demonstrations with behavioral data. The proposed weighting scheme departs from the conventional completeness assumption commonly imposed in the theoretical analysis of deep reinforcement learning. We establish a sharp convergence rate for density ratio estimation and derive the convergence rate for the excess risk of resulting deep $Q^*$ estimator. Extensive empirical evaluations demonstrate that, compared to existing methods, our method achieves significant improvements in numerical performance and policy generalization, providing specific guidance for the rational utilization of expert demonstrations.

Locally stationary Markov chains have attracted a lot of attention in statistics in the past three decades. In this paper we prove a variety of limit theorems for partial sums generated by such chains. We first provide explicit formulas for the asymptotic mean and variance (and also higher moments), and obtain optimal convergence rates towards them. We then prove a Berry--Esseen theorem, a local central limit theorem, Edgeworth expansions, and large and moderate deviations principles. Our approach involves a parametric Perron--Frobenius theorem, which is proved using the theory of complex (Hilbert) projective metrics developed in \cite{Rugh,Dubois}, together with local approximation arguments and ideas in \cite{dolgopyat2023berry}.

Prior information for real-world physical quantities is most elegantly expressed via partial differential equations (PDEs). In this paper, we propose a novel way to solve PDEs using probabilistic inference on a factor graph. In general, factor graphs provide a natural way to encode prior knowledge into a model as explicit factors; here, this knowledge is provided by a governing PDE, which narrows the solution space, while observed data further shape the posterior over the parameters. The approximate parameter posterior is inferred using message passing based on moment matching, without posterior sampling or global gradient-based optimization. We demonstrate our approach on the first-order advection and the second-order semi-linear Fisher-KPP equations, where it achieves predictive accuracy comparable to a standard baseline while providing structured predictive uncertainty. Moreover, the inferred posterior marginal means and uncertainty structure match more closely those obtained using Hamiltonian Monte Carlo than the evaluated mean-field variational inference baseline, while requiring up to 10x less training time in our experiments, with inference speed comparable to variational inference.

A forecasting benchmark reports which method won. We show that the answer is set by the evaluator's choices before any model is fitted. We benchmark 24 forecasting methods and one textbook reference, including six 2025-era time series foundation models, on a production marketplace panel of 1,887 business customers over 67 months. We hold the data, the horizon and the period fixed, and vary only the evaluation design. Three choices each reverse or dissolve a headline conclusion. Changing the unit of analysis from the market total to the individual customer moves our production baseline from second of nineteen, beaten by nothing, to twenty-third of twenty-five. Nineteen of its twenty-four challengers beat it there. Changing how much error is pooled decides whether a Diebold-Mariano test finds anything at all. Scoring prediction intervals rather than point forecasts reorders the field almost completely, with a rank correlation of 0.02 on intermittent demand. We then measure what the deployed system gets from this. Its selection rule captures 55% of the distance between doing nothing and choosing with hindsight. The reversal is not a quirk of our data. We ran the released protocol, unchanged, on the public M5 retail panel. The same baseline shape places first at the market total and last per series, beaten by everything, and a replayed selection rule closes 64.7% of the same floor-to-ceiling distance there. Adding five zero-shot foundation models to that roster changes who wins at the total, not the shape. The bands' blind spot travels too: conformal bands under-cover most on the spikiest items. Splitting our own panel into ever smaller groups turns the contrast into a curve: the baseline's rank worsens at every level of disaggregation. We release the evaluation protocol and report an error of our own that inverted a result before we caught it.

Most changepoint procedures fix segment forms and boundary relationships before estimating changepoint locations. We introduce GRACE, the General Regime-Aware Changepoint Estimator, for an $\ell_0$-penalized problem that jointly selects the number and locations of changepoints, segment forms, transition types, and continuous parameters. A transition depends on the preceding fit only through the quantities it inherits. GRACE retains accumulated cost conditional on these quantities, using joint envelopes when several must be inherited together. This yields an exact conditional-cost calculus for transitions that preserve, release, or fix finitely many boundary parameters. Optimal partitioning, FPOP, and CPOP arise as special cases. CPOP motivates the construction: its level-conditioned cost supports continuous piecewise-linear fitting but cannot represent a discontinuous level shift preserving the incoming slope, because optimality over endpoint level need not imply optimality over slope. GRACE profiles each generated linear candidate at fixed level and at fixed slope before pruning either envelope. The two envelopes support continuous slope changes, level shifts preserving slope, constant regimes, trend termination and resumption, and complete resets within one exact optimization. Under least squares, finite-dimensional linear-basis models with affine boundary relationships yield quadratic candidate costs. We establish exactness of the general recursion and develop functional, penalty-based, and objective-bound pruning rules. Polynomial and seasonal specializations illustrate the framework, whose recursion permits substantial parallel computation.

Diffusion policies represent multimodal action distributions, but an advantage-weighted update does not specify how to sample from the resulting target distribution. We propose Schr\"odinger--F\"ollmer Actor--Critic (SFAC), an offline-to-online reinforcement learning (RL) method for Kullback--Leibler (KL)-regularized policy improvement through conditional diffusion. A minimax Bellman critic estimates the advantage function, which defines an exponentially tilted target policy. A Doob $h$-transform expresses this update as a correction to the reference diffusion drift. We derive a posterior-mean representation of the correction and estimate it using paired self-normalized importance sampling (SNIS). Supervised regression on these drift targets updates the neural actor without critic action gradients. In the small-update regime, the KL-regularized update follows the natural policy-gradient direction, and the Doob correction represents the same local change in the space of diffusion drifts. Under suitable conditions, we derive finite-sample bounds that separate the effects of critic estimation, neural drift regression, finite-sample SNIS, diffusion discretization, and inherited actor error on expected average policy suboptimality. Synthetic experiments assess the accuracy of approximation to prescribed advantage-tilted targets and sensitivity to sampling budgets. On six offline-to-online continuous-control tasks, a reference-anchored implementation achieves higher final-window returns than those of its corresponding offline initialization.

Matrix time series is a series of matrix data observed over time. Analytical tools for such time series is needed in many applications in finance, economics, engineering and many other fields. To avoid the use of vectorization of the matrices which loses the column and row information, and the vector autoregression framework in traditional time series analysis, \cite{chen2021autoregressive} proposed the Matrix Autoregressive (MAR) Model. The model maintains and utilizes the matrix structure, leading to a substantial dimensional reduction and admitting explicit interpretations, comparing with the vector autoregressive model on the vectorized data. However, the MAR model still encounters difficulties in dealing with large dimensional matrix time series as the coefficient matrices in MAR models are also large. In this paper we propose to achieve further dimension reduction through reduced-rank constraints of the coefficient matrices in the MAR model. Estimation and rank determination procedures are studied. Theoretical investigation and empirical examples show that the reduced-rank constraint can achieve higher statistical efficiency than the MAR model.

In causal mediation analysis, natural direct and indirect effects are defined through nested counterfactuals that combine the outcome under one exposure level with the mediator value under another and are therefore inherently cross-world. Their canonical identification additionally relies on cross-world independence assumptions. Consequently, both the estimands and their identifying assumptions remain controversial. In this paper, we ask what additional single-world structure would be required to re-express cross-world estimands as single-world quantities, identifiable under single-world assumptions. We show that these quantities can be written as mixtures of controlled single-world effects under assumptions involving a susceptibility marker, a possibly latent baseline variable that encodes the "would-be" mediator value under a reference treatment arm. If observed, the marker would make several implications of the marker restrictions empirically testable. These results provide a transparent single-world formulation of cross-world mediation effects while making explicit the assumptions required. Although these conditions clarify how cross-world estimands can be interpreted within a single-world framework, their practical relevance depends on whether such susceptibility markers can be justified in specific applications. We discuss implications for principal stratification and show in the Supplementary Material how the formulation extends to ordered mediators and settings with exposure-induced mediator-outcome confounding.

We study online inverse linear optimization with a fixed unknown linear utility: in each round, an environment presents a compact action set, the learner recommends an action from it, and the environment returns an action that maximizes the utility over the same set. When the utility vector and the actions lie in the $d$-dimensional Euclidean unit ball, we give a randomized algorithm whose regret---the cumulative utility shortfall relative to optimal actions---is $O(\sqrt d)$ in expectation for every time horizon, without knowledge of the horizon. The dependence on $d$ is optimal up to a constant factor by the known $\Omega(\sqrt d)$ lower bound for horizons $T\ge d$. Our algorithm maintains matrix multiplicative weights on polynomial feature spaces at geometrically spaced scales. It selects a recommendation distribution by solving a linear program and updates its score matrices by comparing the available actions with the feedback action. With rational oracle outputs and feedback actions, an implementation computable relative to a linear-optimization oracle preserves the $O(\sqrt d)$ regret bound. Whether the same rate is attainable with running time polynomial in the dimension, horizon, and input length remains open.

Vogels et al. (2024) presented an empirical comparison of Bayesian methods for structure learning in undirected Gaussian graphical models. The method implemented in the R package BGGM was run with equal prior probabilities for a null, negative, and positive partial correlation, implying a prior inclusion probability (equivalent to a prior graph density) of 2/3. The other Bayesian methods in the comparison used a prior graph density of 0.2 however. The data-generating densities ranged from 0.01 to 0.1. This resulted in a substantial overestimation of the inclusion probabilities of absent edges by BGGM. For a fair comparison of the performance of the different methods, I repeated the simulation using a prior inclusion probability of 0.2 using BGGM. In this case, the performance of BGGM is comparable with the other Bayesian methods.

Adaptive arm selection changes the distribution of the observations collected by a bandit algorithm, but it need not change their limiting empirical spectrum. We study Gaussian bandit designs in which the dimension and the number of observations grow proportionally. A quantitative coupling theorem compares the design generated by any causal selection rule with an independent Gaussian design. If the logarithm of the number of available arms is sublinear in the dimension, the empirical spectral distribution converges to the Marchenko-Pastur law, uniformly over the selection rule. Consequently, Gaussian Bayesian bandits have policy-independent first-order limits for posterior mean-square uncertainty, squared posterior covariance, and information acquisition. For linear-score selection, we obtain the exact conditional arm distribution and show that two-arm selection produces an exactly Wishart Gram matrix in every dimension, despite its nonzero conditional mean. For a fixed selection direction, we identify an explicit eigenvalue and eigenvector transition governed by the second moment of a Gaussian maximum. A counterexample shows that a direction's overlap with a reference signal does not determine this transition. Finally, an exponentially large arm pool permits a different bulk limit, establishing the order-sharpness of the arm-growth condition. These results distinguish global spectral stability from directional effects in adaptive bandit data.

Researchers increasingly rely on third-party platforms such as Similarweb and Semrush to measure web traffic when first-party analytics are unavailable. Yet these platforms report model-generated estimates rather than raw data, raising questions about whether their measures preserve the temporal and cross-source variation required for causal inference. As a motivating diagnostic, we examine reported referral traffic around two documented search-engine outages; the absence of visible discontinuities illustrates why preservation of identifying variation cannot be taken for granted. We then characterize three mechanisms, within-source smoothing, cross-source leakage, and treatment-induced calibration error, through which platform processing can generate nonclassical outcome measurement error. Analytical results and a stylized difference-in-differences simulation show that this error can attenuate, amplify, or reverse estimated treatment effects. Our findings caution against using third-party traffic measures based on black-box proprietary models for causal inference.

In safe hypothesis testing with test supermartingales, Ville's inequality provides anytime-valid type-I error guarantees for every significance level $\alpha\in(0,1]$, if one rejects the null hypothesis whenever the wealth process first exceeds $1/\alpha$. Due to an inherent asymmetry, the type-II error behaves differently. We prove two things about the latter, for a simple null and alternative. First, the mean growth rate $\mathbb{E}_{P_1}[\log E]$, the e-power, that Kelly betting and growth-rate-optimal e-variables maximise, bounds nothing on its own. For every level $c>0$, every $\alpha$ and horizon $t$ we construct e-variables of conditional e-power exactly $c$ whose probability of not rejecting by $t$ is arbitrarily close to one. It forces eventual rejection, but no finite-horizon guarantee follows. Second, the quantity that does control the type-II error is the Chernoff-Stein exponent of an e-variable, $\Lambda(E)=\sup_{s\ge0}\{-\log \mathbb{E}_{P_1}[E^{-s}]\}$, whose range is exactly determined: $\sup_E \Lambda(E)=\mathrm{KL}(P_0\|P_1)$, the classical Chernoff-Stein exponent, and so the ceiling of its own per-e-variable form. One conditional application of Hoelder's inequality per step gives it, for every test supermartingale on an arbitrary filtered space, with no independence or product structure; the i.i.d. case adds that it is matched, and attained by nothing. The e-power has its own ceiling, $\mathrm{KL}(P_1\|P_0)$, and that one is attained, $P_0$-a.s. uniquely, by the likelihood ratio $R$. The two optima are the same divergence in opposite arguments, at opposite ends of the flattened family $R^{\beta}/\mathbb{E}_{P_0}[R^{\beta}]$: the ceiling as $\beta\downarrow0$, $R$ at $\beta=1$. Which $\beta$ is best is settled by the horizon, exactly: $R$ is optimal at $t=\log(1/\alpha)/\mathrm{KL}(P_1\|P_0)$ alone, beaten by sharpening $(\beta>1)$ below it and by flattening above.

We study structural-change testing in multivariate linear regression when the response dimension is proportional to the sample size and the number of predictors is fixed. Despite its relevance to applications across a broad range of fields, this problem remains underexplored. The alternatives of interest allow multiple unknown changes in a prescribed linear contrast of the predictor effects. We construct a least-squares-based Wald statistic standardized by the residual covariance estimator and scan it over candidate change-point segmentations. We propose flexible segment scans for both single and multiple change points, together with a discretized multiscale scan that reduces the computational cost. When the response dimension and sample size diverge proportionally, we establish weak convergence of the normalized statistic process to a centered Gaussian process, yielding implementable critical values for all three scans. We further characterize their asymptotic power under local alternatives, explicitly describing how the signal, design, and high-dimensional aspect ratio determine the limiting power. The finite-sample performance is examined through simulation studies. We apply the proposed methods to filtered and standardized returns from U.S. equity portfolios to investigate structural changes in their exposures to the Fama--French factors.

Recent literature has shown a strong connection between optimization and sampling. We develop the corresponding first-order theory for diffusion models. First, the SDE-based reverse-time flows of overdamped and underdamped Langevin diffusions contract relative Fisher divergences at explicit exponential rates whenever the stationary potential of the forward process is strongly convex---a condition on the noising process one chooses, not on the data. This is a unique advantage of SDE-based reverse diffusion, absent in the reverse process based on ODEs. Second, we incorporate discretization and establish averaged first-order stationarity bounds---the sampling analog of averaged gradient-norm guarantees in nonconvex optimization---for samplers of both overdamped and underdamped diffusion models. As in nonconvex optimization, the convexity-free certificate is local: it guarantees score consistency, not global mode weights.

This paper introduces a novel semiparametric testing framework for localized point-impact effects within a spatial semi-functional linear regression model. We address an intricate inferential problem that simultaneously accommodates three distinct sources of complexity: an infinite-dimensional global functional slope, an unknown smooth spatial nuisance surface, and data-dependent, unobserved impact locations. To isolate the point-impact coefficients, we construct a robust residualized cross-moment estimator after regularizing the functional component via Tikhonov inverse and estimating the spatial nuisance surface through local-linear smoothing. We establish a multivariate central limit theorem under the null hypothesis of no point-impact effect. We develop two operational test statistics: a general formulation whitened by the long-run covariance matrix to handle arbitrary spatial dependence, and a simpler chi-square calibration valid under short-range spatial dependencies in regression errors. Rigorous finite-sample simulations and an application to Canadian weather data demonstrate the superior power and size control of our methodology.

We study Bayesian linear bandits with an isotropic Gaussian parameter, independent Gaussian candidate arms, and Gaussian reward noise when the horizon is proportional to the dimension. The normalized posterior uncertainty has an explicit limit that is uniform over all causal policies. Gaussian posterior identities then determine the limiting parameter overlaps without an assumed closure of the adaptive recursion. These results yield exact regret curves for Thompson sampling, posterior-mean greedy selection, and a family of policies that scale the posterior sampling covariance. The normalized realized cumulative regret converges in L1, uniformly on compact proportional-time intervals. A policy-uniform lower bound identifies the limiting optimal Bayes regret and proves that posterior-mean greedy selection attains it. Thompson sampling incurs a strictly larger leading regret; its instantaneous regret ratio relative to greedy selection lies between one and two and approaches two at long proportional horizons. Closed-form cumulative curves also identify a different comparison in the vanishing-noise limit. Finally, the instantaneous regret converges to a nondegenerate Gaussian decision-loss distribution, rather than to its mean. The analysis separates the amount of information acquired by a bandit policy from the quality of the decisions made using that information.

This paper studies risk-averse treatment allocation when individuals self-select into treatment based on unobserved characteristics. We develop a framework that combines the marginal treatment effect approach to endogenous selection with a general class of coherent risk measures that capture distributional preferences over welfare outcomes. We show that the planner's problem admits equivalent interpretations in terms of uncertainty aversion, distributional robustness, and worst-case welfare. For law-invariant coherent risk measures, we derive a Kusuoka representation that expresses the planner's objective as a weighted evaluation of different regions of the welfare distribution and characterize the resulting optimal allocation rule. We further establish finite-sample regret guarantees for empirical risk-averse policy learning, showing how the statistical difficulty of learning a policy depends on the planner's sensitivity to adverse welfare outcomes. The framework nests the risk-neutral policy learning model of \cite{Kiatagawa_Tetenov_2018} as a special case. An application to the \cite{Card1995} college proximity data demonstrates that incorporating risk aversion can lead to economically meaningful changes in optimal college admission policies.

Global forecasting models pool many series and learn one shared function. Gradient-boosted trees are their most common form. We measure a failure of this design that has not, to our knowledge, been documented. Train a global tree on the individual series of a hierarchy, then ask it for the hierarchical aggregate. The aggregate sits far outside the model's training range, and the forecast collapses. The model under-predicts the total by 30-50x in our production deployment, and by up to 496x in a public M5 reconstruction. The mechanism is known: beyond its training range, a tree predicts a constant. It surfaces at the aggregate because the total dwarfs every training series. The cure is not new. Per-series scaling, the preprocessing step that Montero-Manso and Hyndman (2021) recommend, prevents the collapse. So do a weighted aggregate-level training row and seasonal differencing. Our contribution is the characterization. The collapse reproduces on five panels: a production business-to-business marketplace, a synthetic hierarchy, M5, Australian Tourism, and a public business-buyer panel. It holds on three tree libraries, is invariant across training seeds, and is statistically significant. Its onset is immediate and tracks a simple support bound: a scale gap of only 1.15x already costs a third of the total. No standard configuration change prevents it: pooling every hierarchy level into training fails at scale, and the one knob that fits linear models in the leaves softens it without curing it. Rolling the forecasts forward recursively separates the cures: the aggregate-row cure re-collapses, per-series scaling degrades but stays low, and only seasonal differencing keeps its one-step accuracy unchanged. We close with a three-step procedure for diagnosing and preventing the failure in deployed systems.

Covariates summarized over a subject's completed follow-up are sometimes entered into Cox regression as though observed at baseline. This practice incorporates future event or censoring information and changes both the estimand and its sampling behavior. We analyze an affine class in which a genuine baseline covariate is contaminated by realized follow-up time. Treating partial likelihood as an observed-data estimation criterion, we derive the population score and characterize its unique least-false coefficient. Exact continuous-time benchmarks show scale reduction and saturation under strong contamination, while administrative censoring destroys the reduction and may produce overshoot. Grouping exit times with Breslow ties changes the geometry: the coefficient has a single hump and eventually returns to zero even though the induced association diverges. We also derive an observed-data influence function and show why model-based variance can be either too small or too large. A subject-level sandwich consistently estimates uncertainty around the least-false target under the stated conditions, but it does not correct the target itself.

Machine Learning and more specifically Deep Learning involves solving large scale nonconvex optimization problems. Several algorithms have been proposed in the literature, that seem to achieve satisfactory practical efficiency for difficult instances, the Stochastic Gradient Method being the most rudimentary, while still outperforming more recent algorithms at a number of learning tasks. A major open question about the current methods used in deep learning is to understand their convergence properties. Following a line of previous works about the long time behavior of gradient-type algorithms, %and in particular the recent contributions from Azizian et al., we present a new approach for studying the asymptotic properties of a wide family of methods from an ergodic theoretical viewpoint. Our main results include a study of the expected time for a stochastic optimisation algorithm to reach a certain small neighborhood of a minimizer and show that this reaching time distributes exponentially around its average, which is given by the inverse of the stationary measure of the target. The assumptions on the Stochastic Gradient noise include the Gaussian and the Sub-Exponential assumptions.

We consider contextual binary prediction with i.i.d. contexts from an unknown distribution and adaptively chosen losses. We show that a simple Follow-the-Perturbed-Leader algorithm with Gaussian perturbation for each observed context achieves the optimal $\widetilde O(\sqrt{T\log N})$ expected regret for a class of $N$ experts, while requiring one optimization-oracle call per round and no explicit enumeration of the class. For an infinite hypothesis class $\mathcal H$, the algorithm attains $\widetilde O(\sqrt{T\operatorname{VC}(\mathcal H)})$ regret. This resolves an open problem posed by Lazaric and Munos (2012), showing that hybrid classification is computationally as easy as statistical learning.

Researchers often summarize mediation in clustered data with a single-level product-of-coefficients estimator. We show what this pooled analysis estimates and why the target depends on the design. The pooled a-path is weighted by the within-between variance composition of the treatment. The pooled b-path uses a different weight, based on the mediator variance that remains after residualizing on treatment and covariates. Because these weights generally differ, the pooled indirect effect can differ from a level-respecting combination of the within- and between-level indirect effects. We derive the pooled probability limit and split this difference into a nonlinearity term and a weight-mismatch term. Equal b-paths remove the distortion, but equal a-paths do not. Allowing the two weights to vary gives a bilinear surface with saddle geometry, sharp corner bounds, and a zero-distortion frontier. We also define a design-standardized indirect effect and practical diagnostics. Simulations verify the identities and show that a common decompose-first, cluster-weighted construction can overstate apparent distortion. In an ECLS-K:2011 illustration with 8,363 children in 753 baseline schools, the within-level weights were 0.723 and 0.896, a difference of -0.173 (95% bootstrap CI [-0.193, -0.151]). The conventional construction produced 1.76 times the absolute distortion of the identity-consistent construction. The identity-consistent AD interval included zero. Thus, the application shows unequal weighting and sensitivity to construction, not evidence of nonzero population distortion.

Online experiments must often be evaluated before long-term outcomes mature. Under rolling enrollment, these outcomes are observed only for early enrollees, while short-term surrogates are available for everyone. We compare seven estimators across eleven data-generating processes, spanning partial mediation, drift, outcome sparsity, and enrollment-time labeling, with up to $R = 2,000$ replications over more than 500 method-by-scenario cells. We find a sharp robustness-efficiency tradeoff: the surrogate index delivers large efficiency gains when surrogacy holds but its coverage collapses under violations, while PPI-family methods stay asymptotically valid under random labeling at smaller gains. We give the finite-sample variance of PPI++ in the all-units parameterization for a fixed predictor, a joint asymptotic distribution for the two estimators under cross-fitting, and a Hausman-type estimator-disagreement diagnostic, then quantify the detection-damage gap: in the partial-mediation design, where we locate both edges, a band of violations destroys surrogate-index coverage yet is too small to detect on most datasets. On the 64,000-customer Hillstrom experiment the diagnostic rarely flags a violation that biases the surrogate index, and a Cauchy-kernel hybrid of the two estimators inherits 10.7% relative bias; on the 14-million-user Criteo experiment the violation is detected, and subsampling traces detection turning on with scale as damage persists. PPI++ has limits: at a rare-conversion $n = 30,000$ Criteo subsample its empirical coverage is 85.0%. We recommend prespecifying PPI++ with the exact variance as the primary analysis under random labeling and adequate labeled outcome counts, reading the diagnostic as a warning, not a certificate, and treating the surrogate index as a sensitivity analysis.

Conventional score-based diffusion models learn scores without representing normalized densities, whereas tractable normalized models support both sampling and direct likelihood evaluation. A recent tensor-network approach provides such a representation but is largely restricted to low-dimensional lattices. Here we introduce a discrete diffusion model that parameterizes normalized probability distributions using variational autoregressive networks. Explicit Markov jump operators govern the forward noising and reverse denoising dynamics, extending discrete diffusion models with normalized distributions to spin systems on higher-dimensional lattices. We apply this framework to the two- and three-dimensional Ising models across ordered, critical, and disordered regimes, accurately computing thermodynamic quantities including free energy, energy, and magnetization. We further integrate the framework with Monte Carlo sampling, using adaptive diffusion steps to maintain high acceptance rates even at low temperatures while enhancing sample diversity. These results establish a neural-network framework for the discrete diffusion model with normalized probability distributions.

Matrix completion, the problem of estimating missing entries in a matrix from noisily observed ones, underlies a diverse array of problems such as recommender systems and counterfactual outcome estimation in panel data. Many algorithms address the problem using regularized least squares, often with the nuclear norm as a regularizer, but this method yields a point estimate with no built-in uncertainty quantification. A Bayesian formulation is a natural alternative, and if the noise variance is known, the nuclear norm-based prior yields a log-concave posterior. Unfortunately, in practice, the noise variance will not be known a priori, so for a fully Bayesian approach, a prior must be imposed on it. We give the first sampler for this model with an explicit non-asymptotic guarantee: polynomial in the matrix dimensions and in the reciprocal of the target accuracy. Our technique is to discretize the distribution of the noise precision onto a grid and build a categorical posterior via thermodynamic integration. This extension is not specific to matrix completion and may be useful in other non-log-concave sampling problems where the non-log-concavity is restricted to a single variable and the joint distribution of the remaining variables is nonsmooth. Our contribution is a feasibility result: we show that a polynomial-time Bayesian sampler for this model exists at all, and the resulting complexity, while polynomial, is not intended as a deployable algorithm at current problem scales.

Variable selection using information criteria has an explicit statistical target but requires discrete search over candidate models. The smooth information criterion (SIC) replaces the discontinuous model-dimension term by a differentiable approximation, with an $\epsilon$-telescoping continuation strategy progressively sharpening this approximation without data-driven selection of a regularisation-strength parameter. We develop SIC as a general procedure for coefficient-level variable selection in generalised linear models (GLMs) and use it to address a central question: how faithfully does smooth optimisation reproduce the corresponding discrete information-criterion selection problem? Focusing on BIC, we benchmark SIC directly against exhaustive subset selection where feasible, using exact support agreement, BIC difference and selection behaviour across a varying signal-strength boundary. Simulations in Gaussian, binomial and Poisson regression show that SIC closely reproduces exhaustive BIC selection and tracks the exact BIC selection boundary. Relative to stepwise BIC, LASSO, SCAD and MCP, SIC produces competitive variable-selection performance while retaining sparse models, with predictive performance broadly comparable across methods. Computational advantages over stepwise BIC increase with predictor dimension. In a real-data application with 16 candidate predictors, SIC recovers the globally BIC-optimal model among all 65,536 supports at a small fraction of the computational cost of exhaustive enumeration.

For a class of regular estimators, H\'{a}jek, in his celebrated ``Convolution Theorem,'' showed that the asymptotic distribution of a regular estimator is the convolution of the distribution of an efficient estimator and some residual distribution. This result constitutes the foundation of the concept of asymptotic efficiency of regular estimators. In this paper, we provide a second-order version of that classical result. Introducing a class of second-order regular estimators with a valid Edgeworth expansion, we derive their asymptotic distribution under contiguous alternatives and show that it is the convolution of the second-order efficient distribution and some second-order residual distribution. This constitutes a second-order extension of H\'{a}jek's convolution theorem. Based on this, we introduce a concept of {\it second-order robustness} for second-order regular estimators. For a class of general Bayes estimators and minimum contrast estimators in time series models, this second-order robustness is used (i) in the characterization of second-order robust priors, (ii) in a comparison between the second-order robustness of maximum likelihood and Whittle estimators.

Time-to-event outcomes are often collected together with auxiliary outcomes that may provide additional information about baseline risk. We develop a semiparametric approach for incorporating such information into Cox regression without specifying a joint likelihood. The survival outcome follows a Cox proportional hazards model and a continuous auxiliary outcome follows a partially linear model, with separate nonparametric baseline covariate effects linked through a quadratic penalty. When these effects coincide and the outcome scores satisfy the information identities and a first-order orthogonality condition, we derive the sandwich covariance of the penalized estimator and show that, for any fixed penalty level, the asymptotic variance of the Cox regression estimator is no greater than that under separate estimation. We further characterize departures from the shared-effect setting. Local differences of order \(n^{-1/2}\) induce an explicit mean shift in the limiting distribution, yielding a direct bias--variance trade-off, whereas fixed differences generally alter the population target under nonvanishing penalization. These results motivate an adaptive penalty that borrows information when the fitted covariate effects are close and approaches separate estimation when a persistent difference is detected. Simulations and a real world data analysis demonstrate the validity and effectiveness of the proposed method.

We develop a model-agnostic framework for noise reduction in high-dimensional time series that explicitly targets optimal recovery of a low-dimensional latent dynamic component contaminated by observational white noise. Under the assumption that the latent dynamics live in a low-dimensional linear dynamic subspace, we characterize the optimal linear projection onto the dynamic subspace and provide a geometric description of the residual error in terms of the relative orientation of the signal and noise spaces. We propose estimators for the dynamic subspace and the optimal projection based on lagged covariance matrices, bootstrap dimension selection, and a low-rank representation of the structured noise. Under mild conditions, the resulting denoised series is shown to converge to its population target at the usual parametric rate. Simulations show that the proposed method can substantially improve subspace estimation, reconstruction error, and one-step-ahead forecast accuracy compared with both orthogonal projection-based denoising and the raw data. The approach is illustrated by empirical applications to high-dimensional stock returns and to a 20-variate time series of macroeconomic indicators.

Non-stationary online learning has attracted much attention in recent years, as static regret is insufficient to guide algorithm design in changing environments. To address this limitation, interval regret and dynamic regret have been introduced as two representative performance metrics that strengthen static regret in complementary directions. Interval regret requires an online algorithm to achieve competitive static regret over every local time interval, whereas dynamic regret evaluates performance against an arbitrary sequence of time-varying comparators. Despite their importance, the relation between these metrics has long remained unclear. Prior work has often regarded interval regret as the stronger notion, based on the intuition that local guarantees should naturally induce global guarantees. Consequently, it is widely conjectured that an algorithm with optimal interval regret should automatically attain optimal dynamic regret. In this paper, we first establish a negative result that refutes this intuition of a metric-level implication. Specifically, for both convex and curved functions (including exp-concave and strongly convex functions), we show that there exist instances in which an algorithm with optimal interval regret nevertheless fails to achieve optimal dynamic regret. We then show how to leverage local adaptivity to obtain optimal dynamic regret. In particular, optimal dynamic regret can be attained by invoking an interval regret minimization process over an enlarged Euclidean ball containing the original convex feasible domain and using a suitable domain-converted surrogate loss. This reduction applies to both convex and curved functions. As a byproduct, we obtain the first proper and efficient algorithm with optimal dynamic regret for exp-concave functions, improving prior results while significantly simplifying the analysis.

We study smoothed online prediction with bounded adversarial responses. This widely studied framework bridges i.i.d. sampling and adversarial covariate selection through a smoothness parameter $\textsf{C}_{\textsf{cov}}$, which bounds conditional covariate densities relative to a fixed, unknown base measure. We propose \textsc{Hedge-Cover}, an information-theoretic algorithm that achieves sublinear regret $\widetilde{O}(\sqrt{\text{Pdim}(\mathcal{F}) \textsf{C}_{\textsf{cov}} T})$ for function classes with bounded pseudo-dimension. The algorithm aggregates a carefully constructed family of experts using \textsc{Hedge}, with a prior that links regret to the number of disagreements between a consistent selector and a target function. We bound this number by exploiting covariate smoothness. This answers an open problem posed in \cite{blanchard2025agnostic} on the minimax optimal adaptive regret of the smoothed online regression problem. We establish a matching lower bound for the class of linear predictors. The main intricacy of the lower bound lies in explicitly constructing a challenging sequential covariate distribution supported on mutually orthogonal hyperplanes. This construction may be of independent technical interest. Finally, we revisit the well-specified setting and quantify the effect of response noise. For conditionally $\nu^2$-subGaussian responses, we extend the existing lower bound under realizable responses by showing that the minimax expected regret is $\Omega((1\vee \nu)\sqrt{(\textsf{C}_{\textsf{cov}}-1)dT})$ for a function class of VC dimension $d$. A corresponding upper bound for ERM matches this dependence on $\nu$.

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

In safe hypothesis testing with test supermartingals, Ville's inequality provides anytime-valid type-I error guarantees for every significance level $\alpha\in(0,1]$, if one rejects the null hypothesis whenever the wealth process first exceeds $\frac{1}{\alpha}$. Due to an inherent asymmetry, the type-II error does not have such guarantees: a heavy concentration of the probability on the lower tail of the log-increments can lead to one catastrophic bet that undoes any amount of accumulated evidence. This paper studies how different hypotheses on those lower-tail probabilities lead to different bounds on the type-II error of the sequential test. They all reduce to one master inequality, which bounds the type-II error at level $\alpha$, at a fixed horizon and sequentially, in terms of a one-sided Legendre transform of the (inverse-)moment generating function of the e-variables, evaluated at one number: the amount by which the lower bound of the accumulated e-powers exceeds $\log\frac{1}{\alpha}$. And, the step is lossless, in the sense, that it extracts exactly a constrained information projection. Every bound presented here is a corollary, obtained by a certain majorant of the above function. The hypotheses are: a finite negative moment; an exponentially small crash probability with a moment on the winning side; a wealth floor with a conditional variance, and its Bernstein variant, which interpolates between a Gaussian regime set by the variance and an exponential one set by the scale; a sub-Gaussian or bounded-tilt lower tail; bounded log-increments; and i.i.d. increments, where the majorant is the truth. We also provide an empirical-Bernstein variant. Each hypothesis may either be read as a condition on the e-variables one has, or as the price of betting with an approximation to the likelihood ratio rather than the ratio itself, which satisfies the weakest condition for free.

Studying heterogeneous treatment effects has become essential in experimental and observational studies. A critical assumption for obtaining reliable treatment effect estimates is overlap, which requires that treated and control units have sufficiently similar covariate distributions. Poor overlap may limit the effectiveness of estimators, especially those based on propensity scores, potentially leading to unreliable results. We investigate the effectiveness of kernel balancing (KBal) (Hazlett, 2020) as an alternative to propensity score methods for conditional average treatment effect (CATE) estimation, particularly in settings with overlap violations. Building on optimization-based balancing approaches, we integrate KBal weights into tree-based methods, specifically, causal forests (Athey et al., 2019) and the X-Learner (XRF) (K\"unzel et al., 2019), to assess their impact on bias reduction and estimation precision. Monte Carlo evidence shows that KBal achieves near-exact balance in a transformed feature space, thereby improving treatment effect estimation in cases where traditional reweighting methods struggle due to extreme weights, finite-sample bias, or insufficient removal of pre-existing confounding bias. We apply the proposed methods to the semi-synthetic IHDP benchmark dataset. Overall, the results indicate that KBal leads to performance improvements, especially in settings with nonlinear treatment effects and limited overlap, making it a useful alternative to propensity score methods.

Evaluating a target policy using data generated by a different behavior policy remains a fundamental challenge in reinforcement learning. While most existing work relies on the standard one-step Bellman residual, we consider a convex combination of one-step and two-step residuals with a fixed mixing weight. In ideal settings, this mixed Bellman formulation can provide a natural bias--variance trade-off between approximation error under a restricted value-function class and the increased variance arising from multi-step importance weighting. To solve this mixed residual optimization, we adopt a minimax formulation involving a critic function. Unlike standard approaches that rely on a fixed functional class, we construct a data-dependent critic representation using predicted future feature directions which effectively induces a kernel adapted to the underlying transition dynamics. This allows the critic to focus on directions that are most relevant for the estimated Bellman error. To control overfitting, we use sample splitting to construct the critic and estimate the value function on separate data subsets. Simulation studies and MetaWorld tasks illustrate the effect of the mixing parameter and show that intermediate residual combinations can improve value estimation in challenging settings.

This paper investigates the numerical resolution of the Beurling-LASSO (BLASSO), a convex optimization framework that promotes sparsity in the space of measures. We consider its application to the estimation of Gaussian mixture models (GMMs) with an unknown number of components and unknown diagonal covariance matrices. Our approach combines the Conic Particle Gradient Descent (CPGD) principle with Riemannian gradient descent, to account for the underlying Fisher-Rao geometry of Gaussian distributions. Our contributions are twofold. First, we provide theoretical guarantees for the convergence of our algorithm. In particular, we establish exponential local convergence under a non-degeneracy condition on the solution and relate this assumption to a separation condition on the underlying statistical target. Second, we address practical implementation aspects of CPGD and present numerical experiments illustrating its performance. On the test cases considered, these experiments suggest that CPGD is more robust to overspecification of the number of components than the EM algorithm. We also investigate the impact of component separation on recovery accuracy.

This article is concerned with the asymptotic shape of quantile surfaces, defined as the set of quantiles at a given level $\alpha$ generated by a controlled one-dimensional distribution. Specifically, when the distribution arises as a linear combination of log-normal random variables and the control is a vector of positive coefficients, we prove that quantile surfaces are globally concave in the left tail ($\alpha\to0$) and globally convex in the right tail ($\alpha\to1$). Moreover, these surfaces exhibit asymptotic separation of scale and shape.

Bipartite A/B tests are experiments in which treatment is randomised over one set of units, while outcomes are measured on another. For example, an online marketplace may test a new pricing algorithm on a random subset of items, while the outcome of interest, say, purchase satisfaction, is measured on customers, each of whom interacts with many items. Existing methods for analysing bipartite experiments assume that every randomisation unit is assigned to treatment or control. In practice, often only a subset participates: platforms cap rollout risk, reserve holdout groups, and split their population across concurrent tests. Ignoring the unassigned units biases estimation, while including them requires care. We construct an unbiased linear estimator for bipartite experiments with partial assignment. Observations must be reweighted not only by the (centred) share of treated connections among participating ones (the exposure) but also by the number of participating connections, each inverse-weighted by its participation propensity, so that units whose connections are well covered carry proportionally more weight. The resulting estimator, EARL (Exposure- and Allocation-Reweighted Linear), is unbiased, consistent, and asymptotically normal; we devise two asymptotic variance estimators and show that it has minimal variance in a natural class of linear estimators. EARL remains unbiased regardless of the experience the unassigned units receive, as long as they contribute to expected outcomes additively; in particular, they can be allocated to other, non-overlapping tests. Our theory is complemented with a simulation study on two public datasets, in which EARL attains up to six times lower error than the strongest existing baseline and, in some configurations, over an order of magnitude lower than Horvitz-Thompson-style alternatives.

Local geometric mixing localizes geometric mixing by requiring geometric convergence to equilibrium in total variation only over finitely many transitions. It accommodates local convergence rates and captures rapid local equilibration, even when global mixing is much slower. We establish and discuss local geometric mixing bounds through Dobrushin contraction. We then apply this approach to Diffusion Path Monte Carlo, a recently proposed Markov chain Monte Carlo method, aimed at leveraging advances in score-based modeling, whose ideal transitions coincide with those of the Proximal Sampler. Our analysis covers both the ideal method and its implementable Metropolis-adjusted counterpart, providing mixing guarantees under minimal assumptions. For the ideal method, these guarantees complement recent spectral gap estimates, which we develop into mixing time bounds.

We study causal discovery where each node is a random function. Previous studies on this topic rely on structural assumptions, e.g., linearity or non-linearity, and distributional assumptions, e.g., Gaussianity or non-Gaussianity. In contrast, we make use of covariance operators to avoid these assumptions. Under functional additive noise models, we propose a new sufficient condition to identify a valid topological ordering based on comparing norms of conditional covariance operators. Taking advantage of this identifiability condition, we develop a new mixed regression model that subsumes linear and non-linear models. Together with variable selection, our procedure yields an estimation of the causal directed acyclic graph (DAG) for functional variables. In theory, we develop the least-squares-type theory of this regression model, and derive asymptotic consistency of order determination, sparse regression, as well as identifying the DAG. Computational algorithms based on discrete observations are provided. Applied to simulated data, our approach performs satisfactorily among existing approaches. A real data example of brain effective connectivity is also presented.

Sequential inference estimates latent states from noisy and incomplete observations. Particle Filters (PFs), a class of Monte Carlo methods based on importance sampling, provide a flexible framework for this task, but often suffer from poor sample efficiency and unfavorable scaling with dimension, partly due to suboptimal proposal distributions. We address these challenges by integrating learned proposals into the PF framework. We introduce Neural Optimal Particle Filters (NOPFs), which learn an amortized approximation to the optimal proposal from offline simulated one-step conditioning tuples. The learned proposal is used as a drop-in replacement in standard PF updates, with samples corrected by standard importance weights so that the method asymptotically targets the same filtering distribution under standard support and density-evaluation assumptions. Across stochastic nonlinear benchmarks of varying inference complexity, NOPFs improve sample efficiency and distributional accuracy over standard PF baselines with modest computational overhead. The approach integrates data-driven proposal learning into classical inference without altering the underlying filtering objective.

Standard diagnostic procedures for assessing the goodness of fit of linear models include tests of the null hypothesis of no residual autocorrelation against the alternative of linear dependence. The literature proposes several portmanteau tests for residual autocorrelation in SARMA models, mainly focusing on two cases: short-term and single-seasonal autocorrelation. Since many time series exhibit multiple seasonal patterns, whose periodic components may interact with one another, diagnostic tests for residual autocorrelation in a multi-seasonal framework are needed. However, the literature still lacks portmanteau tests specifically designed for multiple seasonal autocorrelation. This paper addresses this gap by extending classical portmanteau tests to settings with multiple seasonalities. In addition, the proposed approach jointly tests for both short-term and multiple seasonal residual autocorrelation. Test statistics and their asymptotic distributions are defined. Then, Monte Carlo simulations are used to evaluate the performance of the tests and the consistency with the expected results is assessed using statistical tests. The results suggest that the proposed extension can serve as a useful goodness-of-fit diagnostic tool for the class of mSARIMA models. An application to road traffic data is also provided.

Non-negative matrix factorization (NMF) is a singular statistical model whose Bayesian asymptotics are governed by the real log canonical threshold (RLCT). We study the local geometry of the factorization map and derive an upper bound for the RLCT of NMF. Let $H$ be the model inner dimension and $H_0$ the non-negative rank of the true $M\times N$ matrix. Assuming that the true matrix admits a strictly positive factorization of inner dimension $H_0$ in the interior of the parameter domain, we prove, for smooth positive priors, that $\lambda\leq \{(H-H_0)\min(M,N)+H_0(M+N-H_0)\}/2$. This bound strictly improves the previous bound when $H_0\geq3$. The proof uses a local analytic normal form that separates independent linear coordinates from a residual matrix product. When $H=H_0$ also equals the ordinary rank of the true matrix, we obtain the exact value $\lambda=H_0(M+N-H_0)/2$. Under the standard assumptions of singular learning theory, these results bound the leading coefficients of the expected Bayesian generalization error and the Bayesian free energy.

Chatterjee's estimator uses only two local comparisons per interior observation, leaving a finite-replication variance gap under fixed alternatives. We introduce BRACE, short for blockwise rank aggregation for correlation estimation, which replaces adjacent comparisons with local blocks and averages every within-block rank difference. The block size controls the variance cost of finite local replication. An L2 expansion separates the efficient first-order component of the response-rank statistic from an orthogonal finite-replication component. For an admissible diverging block size, the latter vanishes, so the direct rank estimator attains the information bound. The decomposition yields a consistent variance estimator and Wald confidence intervals under fixed alternatives. At independence, the efficient first-order term vanishes, and the proposed estimator enters a second-order regime in which the same factor controls its null variance.

Neural networks are becoming popular in solving mathematical problems, but stochastic models do not provide mathematical exactness by themselves. This study introduces a discovery-to-certification framework in which high-dimensional variational trial functions are learned in a compact separable representation, spectrally diagnosed and pruned, and then certified exactly through multimodular evaluation. Exact certification makes the numerical proofs fully explicit and independently verifiable. This framework can be run on a standard personal computer. Across three extremal problems, we show that neural optimization can contribute to rigorous mathematics in three distinct ways: by discovering improved constructions, by exposing empirical invariants that lead to proofs, and by revealing optimization barriers whose geometry motivates new analytic or numerical representations. More broadly, these results suggest a path toward AI-assisted mathematics in which flexible computational discovery and exact certification become complementary components of a single rigorous workflow.

Combining $p$-values is a fundamental procedure in global hypothesis testing. In modern high-dimensional settings, however, component $p$-values often exhibit complex dependence and rely on asymptotic approximations rather than exact finite-sample uniform distributions. This paper establishes a unified asymptotic theory for weighted transformation statistics that decouples marginal finite-sample approximation error from the joint dependence structure. We also provide sufficient conditions based on conditional probability bounds to verify the joint-tail conditions. Utilizing this framework, we derive explicit dimension-growth and correlation rates for test statistics operating under asymptotic Gaussian and chi-square calibrations. For non-exact finite-sample statistics, we analyze standardized weighted sums, demonstrating how Cram\'er moderate deviations control relative tail error. Analytical examples demonstrate why both marginal and joint conditions are mathematically indispensable for valid global inference under dependence, and numerical experiments confirm that our asymptotic framework maintains accurate finite-sample size control at extreme significance levels.

Gaussian process stochastic differential equations (GP-SDEs) provide a flexible Bayesian model for unknown continuous-time state dynamics with uncertainty quantification, but learning and inference from noisy and irregular observations remain computationally challenging. To address this issue, we propose GP-SDE Matching, a simulation-free variational framework for Bayesian GP drift learning and continuous-time state smoothing. We analytically marginalize the sparse GP posterior to derive a tractable drift-matching objective that accounts for both the posterior mean and uncertainty of the unknown drift. To handle irregular observations, we further introduce an irregular-time-aware variational state posterior that incorporates the actual observation times during both encoding and continuous-time marginal querying. Experiments on the stochastic Lorenz--63 system demonstrate substantially improved drift recovery and state reconstruction under irregular observations, while five system identification benchmarks show robust forecasting under increasing observation sparsity and competitive performance against existing latent-SDE and state-space methods.

Learning the structure of latent factor models involves two central challenges: (1) estimating the number of latent factors and (2) learning the support of the mapping from latent variables to observed variables. This is especially challenging for nonparametric regimes and nonlinear settings. We propose a method for latent structure learning, based on pairwise dependence measures on the observed variables using a graph-theoretic representation. We show that both the number of latent factors and nonlinear mapping structure can be identified from the distribution of observed variables under mild structural assumptions. Unlike prior work restricted to linear correlations, we establish identifiability and consistency for a general class of dependence measures under nonlinear factor models. This motivates a Dependence Thresholding (DT) algorithm, which jointly estimates the number of latent factors and nonlinear mapping structure from observational data alone. We pair this with a neural network architecture constrained by the nonlinear mapping structure, to recover the nonlinear function. Through simulation studies, we show that the DT algorithm is accurate in practice, even when using flexible methods such as neural networks, and exhibits robustness against violations of its assumptions in high-dimensional settings.

We consider stochastic nonlinear optimization problems with deterministic equality constraints. While unconstrained stochastic optimization is well understood, the interplay between optimality and feasibility in the constrained setting poses significant challenges. Moreover, existing theoretical guarantees for constrained stochastic methods predominantly rely on bounded-variance assumptions, leaving the heavy-tailed noise regime largely unexplored. To address this gap, we propose a novel trust-region method within the stochastic sequential quadratic programming framework, termed TR-SSQP. Our method employs a normal-tangential decomposition in the step computation to balance optimality and feasibility. In addition, we incorporate a normalization mechanism in the design of the trust-region radius, together with Polyak momentum for gradient estimation, ensuring stable updates without gradient clipping. When the trust-region radius and the momentum parameter decay at appropriate rates, we establish global almost-sure convergence of the method. To the best of our knowledge, this is the first asymptotic convergence result for constrained stochastic optimization under heavy-tailed noise. We demonstrate the promising performance of the proposed method through extensive numerical experiments, including comparisons among its variants and with existing constrained stochastic optimization methods.

Self-consistency has become a popular technique for enhancing the reasoning abilities of large language models by generating multiple reasoning paths and selecting the final answer through a majority vote. However, because model providers typically charge users in proportion to the number of reasoning paths generated, they have a financial incentive to artificially increase the path count. In this work, we show that an unfaithful provider can exploit this incentive using a simple, efficient algorithm while avoiding detection by an auditor: by generating and strategically reordering additional reasoning paths, the algorithm makes every path appear necessary to reach the majority. To validate our algorithm, we conduct experiments with multiple instruct models from the Llama and Qwen families, as well as reasoning models distilled from DeepSeek-R1, on benchmark datasets spanning mathematics, science, and question answering. Our results suggest that the distribution of additional reasoning paths generated by our algorithm is heavy-tailed and that substantial capacity to overcharge remains even under the best possible audit designed to keep the false-positive rate below $\alpha = 0.1$.

Large language models are increasingly used as inexpensive judges to evaluate outputs, label data, and assess whether a system meets a desired quality standard. Yet using AI judgments for formal statistical inference is fundamentally different from simply treating them as ground-truth labels: AI evaluations can be biased or noisy, and rigorous hypothesis testing requires explicit control of type-I and type-II errors. We study how to use AI judgments, together with selective human verification, to conduct a valid hypothesis test at minimum cost. We consider a population of items with hidden binary labels. After choosing a fixed pool of items, the decision maker can selectively query AI, send an item directly to a human, escalate an AI-scored item to a human after observing the AI report, or stop once sufficient evidence has accumulated. We derive an information-theoretic lower bound that captures the minimum cost of achieving prescribed testing errors and characterizes the value of AI information and human verification through a report-dependent information frontier. Motivated by this characterization, we develop SCALE, a sequential cost-aware policy that combines selective AI scoring with adaptive human escalation. SCALE is valid at finite sample sizes and matches the lower bound to first order as the target error probabilities vanish. We further extend the framework to an unknown AI-output model using paired AI-human pilot data. Numerically, SCALE approaches Human-only or AI-only testing when one source clearly dominates, while achieving its largest savings when inexpensive AI judgments and selective human verification are both valuable.

An LLM judge scores a bank of responses against a rubric, and the reliability comes back at $0.52$. What has been measured? Judge evaluation has begun borrowing reliability statistics from classical test theory, usually without stating the measurement design each statistic assumes, and we show that three widely portable ones mean something different for a judge than for a test because the judge setting rearranges the roles those designs rest on. First, an internal-consistency coefficient computed over rubric elements contains no scorer facet. Holding one judge's measured error rate fixed at $4.72\%$, KR-20 still ranges from $0.01$ to $0.68$ as the item bank is redesigned around it, and varying judge error moves the coefficient by a comparable amount, so item design and judge error are not separately identified and no single value can be read as a property of the judge. Second, the dependability index $\Phi(\lambda)$ is a ratio of variance components, and the classification probability with which it is sometimes identified differs from it by $0.25$-$0.43$ on our bank and by $0.17$-$0.30$ on simulated data where the underlying model holds exactly. Third, Livingston-Lewis accuracy is indexed to an examinee's own true score on the same instrument, so scoring it against external gold conflates judge unreliability with criterion invalidity. Reviewing the three closest judge-evaluation papers, we found no published instance of these errors, which makes the caution prospective. A coefficient that cannot be attributed to the judge nonetheless travels downstream into deployment decisions and disclosure documents. We therefore close with four reporting lines that keep the attribution attached to the number.

Uniform noise-moment bounds exclude stochastic gradients whose variability increases with the iterate. We study ordinary, single-sample stochastic gradient descent for smooth, lower-bounded, possibly nonconvex objectives under distance-dependent conditional moments. Under second moments alone, a direct descent--displacement argument yields $T^{-1/3}$ expected average squared-gradient stationarity with a horizon-dependent stepsize. An explicit oracle-complexity corollary matches the known smooth Blum--Gladyshev (BG-0) lower bound, including the $Lb_2\Delta^3\varepsilon^{-6}$ and $L\Delta\sigma^2\varepsilon^{-4}$ stochastic terms, where $\Delta$ is the initial objective gap and $\sigma^2+b_2\|x-x_1\|^2$ bounds the variance. Thus unchanged SGD attains the minimax stochastic complexity in this second-moment class. For $p>2$, predictable localization and a Hilbert-space Fuk--Nagaev inequality yield a high-probability bound separating logarithmic variance and polynomial rare-shock contributions. The localization radius is derived from the recursion: no bounded-iterate assumption, clipping, normalization, momentum, or increasing batch size is needed. We also give increasing-confidence rates, an objective-gap-growth refinement recovering root-$T$ stationarity, and stochastic $L^p$-Lipschitz examples. The broad BG-0 optimality statement is distinguished from the smaller mean-square-smooth class, in which additional oracle structure permits faster algorithms.

Tree-Adjoining Grammars (TAGs) have recently been introduced to Nonlinear System Identification (NLSI) as a means of encoding an entire model class as a finite set of grammatical rules, from which candidate models are assembled as trees. Existing TAG-based identifiers rely on evolutionary optimisation and return point estimates of the model structure. This paper instead proposes the TAG framework within a Bayesian setting. A generative prior is defined over tree structures and their parameters, and a Reversible-Jump MCMC sampler with structure-preserving tree moves is used to infer the joint posterior over model structure, parameters and predictions. Two training objectives are considered; that is, a one-step-ahead objective with conjugate parameter proposals, and a simulation-based objective handled by likelihood-free inference. The approach is validated on a simulated polynomial NARX system, the Silverbox benchmark, and wave-loading data from the Christchurch Bay Tower, where embedding Morison's equation as a fixed initial tree yields a grey-box model that outperforms the physics-driven baseline. The results demonstrate that Bayesian TAGs are well suited to quantifying uncertainty in equation discovery for dynamical systems and to fitting physics-informed models.

Modern probabilistic machine learning models increasingly produce multivariate outputs with complex dependence structure, from multi-step time-series forecasts to sample path predictions. Understanding which input features drive the predictive uncertainty is important for risk-aware decisions, model diagnostics, and deciding whether the uncertainty should be mitigated or hedged against. This attribution problem requires a choice of how dependencies between output components are treated. Existing approaches reduce the output to a scalar through aggregation or projection before attribution, thereby obscuring whether features affect marginal uncertainty, dependence structure, or both, while component-wise analyses can miss dependence effects entirely. We close this gap by introducing a hierarchy of three entropy-based Shapley games that make this output-side choice explicit for any ordered multivariate outcome, ranging from per-component marginal entropy to fully joint entropy. The hierarchy isolates a cross-component attribution term that captures how each feature shifts the dependence between output components, a quantity invisible to component-wise methods. We establish a chain-rule decomposition of the joint attribution and characterize the cross-component term through conditional total correlation, providing both closed-form and sample-based estimators. Finally, we demonstrate how the framework captures differences in learned joint structure across probabilistic models from distributional regression to a zero-shot time series foundation model.

Despite the importance of reinforcement learning with verifiable rewards (RLVR), the extent to which it can learn new reasoning capabilities remains debated. Here we study the optimization landscape of RLVR on algorithmic tasks, such as iterated group and quasigroup multiplication. To this end, we map entropy-regularized RLVR over myopic tabular policies onto an energy-based (spin-glass) model over deterministic policies. This mapping upper-bounds what RLVR can achieve, and lets us rigorously characterize the landscape in this tabular setting. We show, both theoretically and experimentally, that for a wide class of models and tasks with uncorrelated inputs, this landscape is benign, containing no local minima that could trap RLVR training. Rather, the practical difficulty of these tasks appears to stem, at least in part, from issues such as diffusive barriers and gradient-estimation error in traversing the landscape. These are genuine obstacles that can prevent a solution from being found, but they are distinct from the landscape itself being rugged. We show that these obstacles can often be mitigated through the choice of entropy regulator. Consistent with this theory, we find that a transformer trained from scratch, using only last-token rewards, successfully learns an algorithmic chain of thought for iterated non-Abelian group multiplications.

Studies enrolling participants on the basis of an extreme baseline value are susceptible to regression to the mean (RTM), such that some observed pre-post change is expected even without treatment. Existing methods estimate and decompose RTM retrospectively. We extend this framework to prospective study design by deriving a closed-form sample-size method for continuous, cutoff-selected, single-arm pre-post studies. The method partitions anticipated total change into that expected from RTM and the residual treatment effect, and powers the study to detect the latter. It uses the general bivariate-normal RTM expression, allowing baseline and follow-up variances to differ, and incorporates the corresponding conditional variance of the pre-post change. To our knowledge, no published method or software implements this cutoff-based RTM framework for prospective closed-form sample-size determination in this setting. The Stata command power onemean_rtm provides solutions for sample size and minimum detectable effect, normal-theory power evaluation, attrition adjustment, and sensitivity analysis. Monte Carlo simulation across 58 scenarios demonstrated accurate Type I error control and showed that achieved power rapidly approached nominal power as sample size increased, with the closed-form calculation generally conservative at very small sample sizes.

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

Formulaic alpha discovery is a core challenge in quantitative trading, as identifying alphas that work well together remains difficult. Recent reinforcement learning (RL) methods formulate this task as a Markov decision process (MDP), but two important issues remain unresolved. First, as the alpha pool evolves, the reward function changes accordingly, making the MDP inherently non-stationary. Second, most existing methods optimize a single objective, typically predictive power, while ignoring other important properties of a high-quality alpha pool. Motivated by these challenges, we propose AlphaPareto, an RL method for formulaic alpha discovery. To address non-stationarity, AlphaPareto augments the state to include both the alpha under construction and the current alpha pool, and applies a large language model (LLM) to encode the pool. This design allows the agent to adapt to the evolving search environment. To overcome the limitation of single-objective reward design, AlphaPareto replaces the scalar reward with a multi-objective vector-valued reward that simultaneously captures predictive power, temporal stability, perturbation robustness, and diversity, and optimizes these objectives through a Pareto-regularized learning procedure. Empirical applications to real-world datasets show that our AlphaPareto method outperforms its competitors.

We present spsurv, an R package for semi-parametric time-to-event regression based on Bernstein-polynomial estimation of unknown baseline functions. The package provides a unified modelling interface for proportional hazards (PH), proportional odds (PO), and accelerated failure time (AFT) models for right-censored data, with either maximum likelihood or Bayesian estimation via Stan. Smooth baseline hazard, odds-function, or log-time structures are estimated without assuming a parametric baseline family, while retaining familiar hazard-ratio, odds-ratio, and time-ratio interpretations. We describe methodology, implementation, and syntax; evaluate finite-sample behaviour in a Monte Carlo study; and illustrate usage with oncology trials.

Validating a model on a time series asks for three things at once: each training run should use most of the sample (sufficiency), the test sets should together cover most of the sample (coverage), and training data should come before test data (causality). We prove that the three cannot be had together and price each one. Let $\alpha$ be the smallest training fraction over folds, $\beta$ the fraction of the sample covered by tests, $\Lambda$ the fraction of the sample used as training data from the future of a test point, and $\delta$ the distance from a test point to the nearest training point in its future. Every scheme on a sample of length $T$ satisfies $\alpha+\beta \le 1+\Lambda$ and $\alpha+\min\{\beta,\delta/T\} \le 1$, and under $\beta$-mixing the leakage bias at a test point is at most $2M\beta_{\mathrm{mix}}(\delta)$. In words: going beyond the causal frontier $\alpha+\beta=1$ requires training on the future; that future data must sit within $(1-\alpha)T$ of a test point; and its harm depends on its distance, not its amount. Hence expanding walk-forward is exactly the Pareto frontier of causal validation, $k$-fold cross-validation buys the most future data, and purged $k$-fold with an embargo pays in distance instead, which is cheap when the process forgets quickly but cannot repair the part of causality demanded by non-stationarity. On pure noise, shuffled 5-fold reports an information coefficient of $+0.32$, while contiguous 5-fold, using the same amount of future data, reports $+0.004$.

We introduce Rolling Conformal Prediction (rolling-CP), a distribution-free predictive inference method for the setting of sequential model training. Specifically, given a data stream $(X_1,Y_1),(X_2,Y_2),\dots$, at each time $n$ the trained model may depend on the observed history $\{(X_i,Y_i)\}_{i<n}$. This setting arises naturally in modern sequential training, including one-pass training over massive datasets and continual fine-tuning or test-time adaptation of language models during deployment. Rolling-CP first calibrates each incoming observation against the current predictor and then rolls it into future training. In this way, we avoid the need for data splitting. Remarkably, although the models at times $n=1,2,\dots$ may have entirely different properties and accuracy levels, for exchangeable data it is nonetheless possible to establish a guarantee of marginal coverage, with a familiar universal factor-two guarantee (a worst case guarantee of $1-2\alpha$ coverage, as compared to the target level $1-\alpha$), without any assumptions of stability or any restrictions on the model training process. For i.i.d. data streams, we further prove high-probability training-conditional validity uniformly over time; under stability conditions, coverage guarantees sharpen towards $1-\alpha$. Numerical experiments on sequential regression, multiclass SGD, and one-pass neural-network training further demonstrate the practical effectiveness of rolling-CP.

Objectives: There has been increased appreciation of the differences between marginal and conditional estimates and different types of effect measures regarding their applicability to different target populations. This issue of transportability is of concern in model-based cost-effectiveness analysis (CEA) when treatment effects from (international) trials are applied to (country-specific) baseline risk estimates. The objective of this paper is to create awareness regarding the issues that arise when using different types of treatment effect and baseline risk estimates in a model-based CEA to inform health technology assessment (HTA). Methods: We clarify collapsibility, marginal versus conditional estimation, and transportability; derive the ideal modeling approach implied by a marginal cost-effectiveness estimand; and characterize the issues of common modeling approaches, illustrated with a fictitious state-transition model. Results: An individual-level simulation that predicts outcomes from conditional inputs and averages them over the target population targets the marginal cost-effectiveness estimand. Cohort-model approaches that marginalize inputs early, evaluate an outcome regression model at mean covariates, or combine a conditional effect with a marginal baseline (or vice versa) can misstate cost-effectiveness results even with correct-population inputs; inputs from the wrong population add further error. Conclusions: The most rigorous approach is an individual-level simulation that carries conditional inputs (baseline risk, treatment effect, prognostic effects and effect modifiers) and marginalizes late. Cohort approaches instead marginalize early, relying on aggregated inputs, and do not necessarily target the marginal cost-effectiveness estimand of interest for HTA. Model developers should document, for each input, whether it is marginal or conditional and its population.

Covariates are widely used in randomized experiments to improve precision. However, in the presence of interference, where outcomes may depend on the treatment assignments of other units, standard covariate adjustment methods may fail to preserve desirable properties such as the no-harm property, meaning that incorporating covariates does not worsen estimator performance. Existing approaches that incorporate covariates under interference typically rely on specifying a particular interference structure, and their guarantees can be sensitive to misspecification. In this paper, we study how to incorporate covariate information in a way that preserves the no-harm property while remaining largely agnostic to the underlying interference structure. We focus on the estimation of the expected average treatment effect (EATE) using the H\'ajek estimator and study rerandomization under interference, a design-stage procedure that restricts the assignment space to allocations with sufficiently small covariate imbalance. Our framework also allows the covariates used for rerandomization to depend on the treatment assignment itself, such as the proportion of treated neighbors, which naturally arises in settings with interference. We show that, even under interference, rerandomization can improve estimation precision asymptotically relative to unrestricted Bernoulli randomization, relying only on mild conditions on the dependence structure across units. When a conservative dependence graph is available, we further develop an optimization-based conservative variance estimator for inference under rerandomization.

Metric graphs are generalisations of linear networks and provide a natural framework for the definition of continuously-indexed Gaussian processes. We define a new class of distances on these topologies, termed polyharmonic distances, which unify and extend the spectral construction underlying the effective resistance distance and the biharmonic one. We give both a spectral and a variational characterisation. Furthermore, we show an explicit class of stochastic processes whose variograms coincide with the squared polyharmonic distances. Finally, we show how these metrics can be composed with suitable completely monotonic functions to define isotropic processes having any prescribed finite-order mean-square differentiability along the edges and satisfying the Kirchhoff conditions up to order one at the vertices.

We study conformalized quantile regression and conformalized median regression that fit a model on one block of a time series and calibrate the conformal interval on the adjacent block. The existing theory of conformal prediction for time series rests largely on mixing conditions, which are hard to verify from a time-series model and fail for many standard processes, including simple ones with short memory. We replace this theoretical toolbox with the functional dependence measure, which in principle accommodates long-memory observations. The accuracy of the conformal interval length for time series has been understudied. To the best of our knowledge, this paper is the first work that establishes non-asymptotic coverage guarantees and accuracy of interval length simultaneously for split conformal regression on time series. Furthermore, for Gaussian linear processes with long memory, where both the estimation of the center and its calibration converge slowly, we establish a sharper rate for the length error. We show that the calibrated length converges faster than the estimated center itself, and provide a matching lower bound for the usual centers when the calibration block is sufficiently large relative to the training block. To our knowledge, this is the first theoretical analysis of conformal interval length dedicated to long memory.

We propose a new $L_2$-type test for white noise which allows the dimension $p$ of the time series to either (i) be a fixed constant, or (ii) diverge with the sample size $n$. The proposed test statistic exhibits an interesting phase transition, following two different regimes of behavior: $p$ is fixed, and $p\rightarrow\infty$. Because identification of the operable regime is difficult, if not impossible in practice, we devise a novel adaptive bootstrap method to construct unified testing procedure across different phases. Numerical experiments confirm the good finite sample performance of the proposed adaptive $L_2$-type test in comparison to the existing methods in the literature. The proposed testing procedure has been implemented in R package HDTSA.

We study sequential testing of a fixed family of hypotheses when observations are costly and a sampling horizon is specified in advance. The challenge is to pool evidence for earlier decisions when the number and identities of false hypotheses are unknown, while controlling the probability of any false rejection at level $\alpha$. Existing merges attain the pooled growth rate at a single number of false hypotheses: averaging when one is false, multiplying when all are. Under an independent-stream model with common simple null and alternative distributions, we test each intersection with a prior-weighted mixture of products of marginal likelihood ratios. Closed testing combines these elementary-symmetric-polynomial mixtures to identify individual false hypotheses. Design-specific boundary calibration gives finite-horizon family-wise error control, with exact finite-state guarantees or a confidence qualification for Monte Carlo calibration. The prior-matched mixture uniquely maximizes expected log evidence at each horizon. Mixtures assigning positive weight to every nonempty subset of streams attain log-growth rate $lD$ when $l$ streams follow the alternative. Here $D$ is the mean log likelihood ratio per alternative observation, and a round supplies one observation per stream. This rate attains the first-order intersection-delay lower bound as $\alpha\downarrow0$ at fixed dimension, configuration, weights, and a long enough horizon. Power for an individual hypothesis cannot exceed the best single-stream power at a given deadline, but closure removes the multiplicity penalty when all are false. Gaussian, basket-trial, language-model and advertising studies illustrate both. The primary basket boundaries are 40-52% below $1/\alpha$. Across 41 simulated configurations, the median reduction in capped mean patient outcomes relative to prespecified interim-look Bonferroni tests is 31%.

Conformal prediction provides distribution-free finite-sample marginal coverage, but post-hoc calibration data may be too scarce to learn how uncertainty varies across inputs. Meanwhile, abundant covariates can often be labeled cheaply by domain models or general-purpose language models. We study whether these synthetic labels can improve conditional coverage when only a small trusted sample is available. Building on score-quantile regression, we introduce prediction-powered quantile learning: a synthetic-labeled pool estimates pinball risk, paired trusted and synthetic outcomes correct its bias, and an independent trusted split performs final conformalization. Profiling pinball risk over scalar corrections reveals that population conditional-coverage error is its functional gradient; the corresponding Hessian removes global shifts and weights remaining shape error by boundary density. Composing this geometry with prediction-powered learning yields a three-resource expansion and a benefit--cost rule for synthetic power. Across eight regression benchmarks, synthetic-powered quantile learning substantially improves downstream conditional coverage while preserving marginal validity and producing more compact prediction sets. A human-rating study finds similar gains from external LLM labels and exposes a quality--quantity--cost tradeoff.

We study asymmetric simple binary hypothesis testing between $H_0:P_0^{n}$ and $H_1:P_1^{n}$, based on $n$ independent and identically distributed observations. Leveraging a variational representation of R\'enyi divergence of order $\alpha$, we derive our main result: a finite-sample converse with $\alpha>1$. The bound uses both directions of the divergence $D_\alpha(P_1\|P_0)$ and $D_\alpha(P_0\|P_1)$, tensorises under product measures, and contains familiar data-processing converses as boundary cases. For comparison, we apply the same variational approach to general $f$-divergences and specialise it to total variation, $E_\gamma$, Hellinger, and Kullback Leibler divergences, thereby recovering familiar converses within a unified framework. Together with an achievability bound involving R\'enyi divergence with $\alpha\in (0,1)$, the main converse recovers the phase transition of the optimal Type II error under the exponentially decaying Type I error constraint $\varepsilon_n=e^{-nr}$. Under regularity conditions, the optimal Type II error vanishes exponentially when $r<D(P_1\|P_0)$ and converges exponentially fast to one when $r>D(P_1\|P_0)$. We also derive sample-complexity bounds and extend both the converse and achievability analyses to locally differentially private observations, quantifying the cost of privacy and recovering the non-private achievability bound as the privacy constraint vanishes.

We propose a class of extended generalized Pareto models for right-censored survival data, with particular emphasis on tail inference and long-term extrapolation. Our framework integrates extreme value theory and survival analysis, combining a generalized Pareto distribution with a flexible perturbation distribution on the unit interval. We consider three perturbation specifications: a parametric Beta model, a Bernstein polynomial estimator, and a histogram-based estimator. To facilitate direct comparison, all three models are fitted using a unified iterative procedure adapted to right censoring through Kaplan-Meier-based pseudo-observations. A Monte Carlo study evaluates finite-sample performance across different tail indices, censoring levels, sample sizes, and model complexities. The results reveal a trade-off between flexibility and stability: the Beta specification generally performs best for tail-index estimation, whereas the histogram estimator performs particularly well for scale estimation under low censoring. The Bernstein estimator shows intermediate performance and greater sensitivity to sample size and censoring. Applications to bladder cancer recurrence and heart-failure survival data show that models with very similar in-sample fits can nevertheless produce markedly different tail-index estimates and long-term extrapolations. These findings emphasize the importance of perturbation specification when extended generalized Pareto models are used for survival extrapolation under censoring.

Worst-case Optimal Recovery (OR) and Bayesian learning describe the same Gaussian-quadratic-Hilbert problems in two vocabularies. We sharpen the correspondence - the radius of information equals a nugget-optimized GP posterior variance and is attained by the posterior mean at a closed-form balance nugget - and measure, inside three published Bayesian systems, where the worst-case side pays. The ledger is two-sided: the losses instruct as much as the wins. Morozov calibration tracks a test-access oracle within $1.00$-$1.19\times$ where $\sigma$-blind rules fail, is $4.9$-$6.3\times$ more reproducible across noise draws ($p=0.002$-$0.004$), and is the only deployable rule whose selection survives a change of backend ($1.36\times$ against $12$-$30\times$ for the released weight, ML-II and GCV); tight certificates cover at the information-theoretic floor with no numerical slack. But on exchangeable data split-conformal beats the OR head on interval score, a water-filling prior adds nothing without an oracle noise hint, and under covariate shift the OR band keeps coverage on every dataset yet loses interval score to split-conformal, and to a feature-free constant band, on most cells; what pays is not shift but shift on a learnable target, which a training-free audit statistic predicts before any model is fitted. In Bayesian optimization the certified width is a validity floor whose scalar inflation we prove inert under a checkable margin condition and check at every step. Inertness is graded, not binary, and in the size of the inflation as much as in the objective: $\kappa{=}2$ is inert wherever $\kappa{=}5$ is and on more cells besides, while $\kappa{=}5$ moves half the Ackley seeds and every Griewank seed. Exploration is a shape problem, not a scale one. The design rule: match the guarantee tool to the data regime, and audit the regime first.

Under Gaussian sampling, removing location and scale turns the sample into a direction that is uniformly distributed on a sphere of dimension $n-2$, and sample skewness and kurtosis become a cubic-quartic polynomial image of spherical measure. We use this geometry to obtain exact finite-sample results for the Jarque-Bera statistic $\mathrm{JB}_n$. First, for every $n\ge 3$ the largest attainable value of $\mathrm{JB}_n$ is $n\{(n-2)^4+4(n-1)^2\}/\{24(n-1)^2\}$, attained by a single-outlier sample, so the asymptotic $10\%$ and $5\%$ tests have size zero for $n\le 6$ and the $1\%$ test for $n\le 7$; near this maximum the density of $\mathrm{JB}_n$ behaves like an explicit multiple of $(J_n^+-x)^{(n-4)/2}$. Second, passing from residual coordinates to power sums writes the joint density of skewness and kurtosis as an integral of the reciprocal square root of a polynomial discriminant; for every $n\ge 5$ the innermost integral runs over a single interval and is a Lauricella $F_D$ period. Third, this yields explicit laws: an arcsine law for $n=3$, an algebraic joint density for $n=4$, and a single complete elliptic integral, equivalently ${}_2F_1(\tfrac12,\tfrac12;1;\cdot)$, for $n=5$. Residual-coordinate collisions generate the discriminant singularities of the joint law, while stationary points of $\mathrm{JB}_n$ on the residual sphere govern the singularities of its one-dimensional density.

This note presents a general form for a proper scoring rule that can be used for possibly right-censored survival data without requiring any knowledge or assumptions about the censoring process.

This paper develops prediction models for population-level survival curves observed over time and sampled from a heterogeneous mix of populations. We consider a discrete-time setting where each curve is only partially observed and forecasts of the remaining trajectory are needed for downstream decision making. Our approach recasts cross-population heterogeneity into a multivariate sampling model. We propose two forecasting models for partially observed curves: a full factor analysis model that extends a general factor representation to incorporate the partially observed survival curve, and a double PCA model. The methodology is motivated by demand forecasting in short-term rental markets, where market-level occupancy paths can be viewed as survival curves over the booking horizon and where forecasts of future occupancy feed into dynamic pricing algorithms. We apply the models to the newly released Wheelhouse dataset, which contains time series of market occupancy curves for 500 markets from 2017 to 2022. Model performance is assessed using the integrated quadratic distance, and we compare the proposed PCA-based methods to Holt's linear trend model across multiple forecast horizons. The results show that the proposed models yield accurate and stable forecasts of the remaining survival trajectory and generally outperform Holt's method, particularly at longer horizons.

Large language models (LLMs) as a judge enable scalable evaluation, but their judgments can be sensitive to response order and, even after removing such position effects, can still diverge systematically from human preferences.We introduce DIAL, a unified framework that combines abundant LLM comparisons with limited human comparisons to separate judge-specific position effects, learn shared structure in position-debiased LLM preferences, and adaptively calibrate that structure toward the human preference target. Theoretically, we study three aspects of DIAL: (i) identification of latent LLM preferences, position effects, and human calibration; (ii) adaptive estimation that balances LLM anchoring against limited human evidence; and (iii) fixed-weight uncertainty quantification for the calibrated human preference. Empirically, we evaluate position debiasing and human alignment separately in controlled simulations and on three human-preference benchmarks, showing that DIAL remains robust to unbalanced response order, achieves strong human-aligned rankings with limited labels, and adapts toward human evidence when LLM information is imperfect. Our real-data study collects over 410K judgments from 21 LLM judges in both display orders, providing a resource for future studies of LLM-judge bias, heterogeneity, and human alignment.

We investigate $c$-interior pointwise maximal leakage (PML) as a tool for contraction analyses and disclosure control. Based on the strong adversarial threat models from maximal leakage, $c$-interior PML generalizes local differential privacy (LDP) to data-generating distributions with densities uniformly bounded away from zero by $c>0$. Viewing $c$-interior PML as an algebraic constraint on a kernel yields more flexible (and often tighter) contraction analyses than standard LDP. We provide tight bounds on the Dobrushin coefficient, and bound the contraction coefficient of the Hockeystick-divergence. We further derive strong data processing inequalities on $f$-divergences under $c$-interior PML constraints when the input distributions to the divergence are restricted to be in the $c$-interior. These results extend beyond the regime of pure LDP to cover a larger class of kernels, including, e.g., arbitrary stochastic matrices. We apply the results to minimax theory and provide asymptotically optimal strategies under $c$-interior PML constraints for binary hypothesis testing and mean estimation. The results show that disclosure control with PML allows analysts to reason about systems in a more differentiated manner: For example, it allows us to quantify the privacy leakage of deterministic systems, and can give precise adversarial guarantees with respect to arbitrary distributional assumptions. Interestingly, a recurring theme in the disclosure analyses is that if the privacy problem is relatively regular (if the density bound $c$ is large), private inference can be possible without incurring any additional cost in terms of sample complexity.

We initiate the study of approximating the top eigenvalue and eigenvector of a random symmetric matrix $ A \in \mathbb{R}^{n\times n} $ using $ q(A)b $ where $q$ is a degree-$d$ polynomial and $b$ is a standard Gaussian vector independent of $A$. For spiked GOE $ Y = \lambda vv^\top + X $, we identify $ d_\star = \frac{\log(n)}{2\log(\lambda)} $ to be the critical degree threshold above which accurate approximation of the top eigenvalue and eigenvector is possible. This sharpens the common belief that spectral methods can be implemented by $ O(\log(n)) $-step power iterations and offers a precise connection between spectral methods and low-degree polynomial algorithms, a popular proxy for all polynomial-time algorithms. For GOE $X$, we identify $ d_\star = n^{1/3+o(1)} $ to be the critical degree threshold for top eigenvector approximation, whereas constant degree suffices for top eigenvalue approximation. Moreover, in the limit where $ d/n^{1/3} $ converges to a positive finite constant, we compute the exact asymptotic eigenvector approximation accuracy in terms of the expected squared overlap. These results significantly improve upon predictions made in randomized numerical linear algebra for deterministic data matrices that the iteration count of power methods is governed by the inverse spectral gap. Technically, our analyses leverage extremal properties of Chebyshev polynomials and draw upon the rich literature of random matrix theory.

We study the impact of conditional complier average causal effect (CCACE) estimation methods on the performance of subgroup discovery and heterogeneous causal effect estimation under imperfect compliance. Building on the Bayesian Causal Forest with Instrumental Variable (BCF-IV) (Bargagli-Stoffi et al. (2022)) method, we introduce a two-step, model-agnostic approach that allows any suitable machine learning method to be used for the CCACE estimation in the first step. Specifically, we implement two non-Bayesian tree-based methods, both using forest-based learners: DRRF-IV, a debiased transformed-outcome regression-forest approach, and a GRF-based IV adaptation of the generalized random forest framework (Athey et al., 2019). Through a simulation study, we assess the precision, bias, and the ability to correctly identify the underlying subgroup structure of the proposed methods relative to BCF-IV. The results show that the non-Bayesian methods perform competitively across the considered simulation settings, with performance improving for larger sample sizes and moderately large treatment effects, while reducing computational runtime. We apply our new methods by revisiting an empirical study that examines the effect of prompt admission to intensive care units (ICU) on 28-day mortality across 48 UK National Health Service hospitals. While previous work finds no significant overall treatment effect, we investigate whether subgroups of patients may benefit more from prompt ICU admission.

To design systems that are protected against events much rarer than the observational record, extreme-value methods are needed to extrapolate distribution tails. Tail-index estimators such as the Hill estimator are central to this extrapolation, but overestimating the tail exponent can lead to substantial underestimation of rare-event probabilities. Motivated by this, we derive large- and moderate-deviation asymptotics for the Hill estimator and use them to construct estimators whose probability of exceeding the true tail index decays at a controlled exponential rate (the decay rate). In the large-deviations regime, we show that a simple rescaled version of the Hill estimator achieves an optimal balance between bias and decay rate among scale-invariant estimators based on the same top $k$ order statistics. Under a second-order condition, we quantify the effect of the Hill bias, analyze a bias-corrected estimator, and identify sufficient conditions for moderate deviations in the boundary case where the second-order parameter $\rho$ equals zero.

Leaderboards rank models by their average scores on benchmark items, and they are consulted repeatedly while the evaluation is still running. Existing confidence intervals for a model's rank control their error rate only if they are computed once, after a number of items chosen in advance. If they are recomputed as results arrive, and the evaluation stops once they look decisive, their error rate exceeds its nominal level. Anytime-valid methods keep their guarantees at all sample sizes simultaneously and hence under any stopping rule. They exist for the accuracy of one model, for one pair of models and for the set of models that may be best. For pairwise battles they also give ranks. None gives ranks when all models are scored on the same items, which makes their scores dependent. We construct rank confidence sequences: for every model, a set of ranks that contains its true rank, simultaneously for all models and at all times, at a chosen error level $\alpha$, in finite samples. The construction combines betting e-processes, one for each ordered pair of models, with closed testing over the possible orderings of the models. It allows any dependence between the models' scores on an item. The method has two advantages. A leaderboard can be inspected after every item without inflating its error rate. The evaluation of each model can stop as soon as the question asked about it is answered, which saves compute. When results are examined only once, halfway through or later, little power is lost relative to fixed-sample methods. The paper quantifies these advantages in simulations and on public leaderboard data.

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

Classical continuous-time non-linear filtering fails in generic non-linear state spaces due to an infinite Lie algebraic derivative explosion ($\dim(\mathcal{E}) = \infty$), leading to filter divergence. To resolve this four-decade crisis, we introduce the SMG-Yau-Yau filter and lossless data assimilation framework built on Statistical Fiber Theory over an Orlicz manifold $\mathcal{M}$. By equipping $\mathcal{M}$ with a Riemannian submersion and an Ehresmann connection, the unconstrained Duncan-Mortensen-Zakai score velocity field is orthogonally decomposed into Statistically Verifiable Directions ($\text{SVD}\chi_f$) and Structural Internal Directions ($\text{SID}_f$). System non-linearities and unclosed Lie commutators are orthogonally quarantined in $\text{SID}_f$, protecting macroscopic base parameters from spatial derivative pollution while preserving total score variance energy. We unify the asymptotics through a Dual-Axis Collapse Mechanism: proving our model identically recovers classical Yau-Yau dynamics when $\dim(\mathcal{E}) < \infty$, while large-sample limits ($N \to \infty$) induce a thermodynamic quench that flattens infinite-dimensional geometry under $\dim(\mathcal{E}) = \infty$. Finally, we formulate the Active Acausal Tension (AAT) functional to monitor accumulated model misspecification online. Exceeding a topological capacity threshold triggers Gauge Symmetry Breaking (GSB), which dynamically expands base coordinates ($d \to d+1$) to ensure non-asymptotic stability and convergence to the exact state density.

The analysis of experimental or observational data often requires controlling for the effects of variables of secondary interest. Two common approaches to account for such parameters are grouping experimental or observational units into homogeneous groups (blocking), and observing quantitative variables on each one of them (covariance analysis). Often, both can be used and; therefore, comparing their performance poses a research question. Specifically, under randomized complete block designs where the blocks are defined by quantitative variables, there is interest in comparing the two-way Anova (TWA) with analysis of covariance (Ancova) models. Motivated by real-life problems coming from statistical consulting, we compared the performance of linear models considering block effects and those considering covariate effects from analytical and empirical perspectives. Models assumed Gaussian errors and included the effects of T treatments. The formal comparison relied on the theory of orthogonal projections and vector spaces under containment relationships of column spaces of the design matrices induced by conditions on the experimental design or data structure. We found asymptotic and small-sample conditions for the TWA to be at least as precise as the Ancova, as well as for these models to be equivalent. As to the empirical perspective, models were compared through a simulation study motivated by agronomical experiments. The results suggested that when the covariates show complex heterogeneity patterns, the Ancova outperforms the TWA with blocks defined by traditional approaches in terms of precision of location parameters, but not in hypothesis testing accuracy.

In this paper, one provides a comprehensive mathematical and practical overview of Functional Data Analysis (FDA) specifically applied to parametrized curves in Rp. One observes that curves depending on continuously parameter are naturally present across many fields, such as, for example, monitoring child development in pediatrics, assessing meteorological phenomena, analyzing financial portfolios, tracking neurological functions, and mapping geographic pollution levels. One considers a formal theoretical framework given by a Cartesian product of real separable Hilbert spaces, to model these curves mathematically. To bridge the gap between discrete experimental measurements subject to errors, and smooth continuous functions, the paper details the essential phase of data smoothing and fitting and sets its mathematical formulation via of the Ordinary or Penalized Least Squares criteria which, for the later, using the Sobolev spaces framework, incorporates a smoothing parameter ${\lambda}$ and differential operators to prevent erratic geometric behaviors by penalizing excessive curve roughness. Furthermore, one extends the usual Principal Component Analysis to its functional counterpart in $\mathbb{R}^3$, using the calculus of variations and the Euler-Lagrange multiplier theorem in order to find eigenfunctions and eigenvalues of the covariance operator to exhibit the optimal decomposition of the spatial variance and finally, one highlights the geometric advantages of FDA over traditional multivariate data analysis, emphasizing its unique capacity to capture critical differential features as velocity, acceleration, curvature, and torsion.

Language model post-training is often bottlenecked by the need for human-collected preference data, which is expensive and difficult to scale. Reinforcement learning from AI feedback (RLAIF) style approaches that leverage pseudo labels offer an abundant alternative but introduce systematic biases that degrade downstream alignment. Recent general-purpose semi-supervised methods correct for teacher bias using a small set of human-labeled examples, but suffer from high variance especially when human annotations are scarce. To this end, we propose ABC-Align, leveraging abundant pseudo label signal to minimize variance and applying a lightweight, adaptive correction grounded in the human-labeled subset. The correction strength is tuned automatically during training using plug-in estimates of the relevant bias--variance quantities. On LLM alignment with RLHF, DPO, and GRPO where human feedback is scarce, we empirically demonstrate that ABC-Align achieves superior performance over prior semi-supervised baselines in a series of experiments on an increasing scale. Our code is available at https://github.com/SewoongLab/abc-align .

Neural Processes (NPs) are model-based meta-learners that implicitly learn a stochastic process and adapt to a new task from a small context set. Most extensions of NPs focus on improving the neural network architecture. We instead develop an extension motivated by the shared hierarchical interpretation of meta-learning and mixed-effects models. Specifically, we introduce Neural Process Boosting (NPBoost), which decomposes structured response variability into tree-boosted fixed effects shared across tasks and NP random effects that capture stochastic task-to-task variation. We propose to train the two components jointly using a boosting algorithm in which an NP learns residual task-specific structure and a tree ensemble estimates common patterns across tasks. Across synthetic and real-world tabular meta-learning problems, this decomposition improves over a standard NP when the shared structure contains discontinuities or other irregular patterns that boosted trees can represent effectively.

In target trial emulation, the treatment strategies that patients follow are inferred from the treatments they actually receive in routine care. However, in most settings, the outcome may preclude the observation of planned treatment, giving rise to immortal time bias through misclassification of treatment strategy, while markers of treatment response may influence subsequent treatment decisions, giving rise to time-varying confounding. A key step toward unbiased treatment effect estimation is to partition follow-up into sufficiently short time intervals to unfold the feedback relationships involving treatment and represent the resulting causal relations with a directed acyclic graph. In this study, we present the possible within-interval causal orderings induced by this partitioning, discuss their causal implications, and assess their plausibility across clinical settings. For each causal ordering, we derive the corresponding g-formula. Using ancestral multi-world networks and simulations, we show that the standard cloning-censoring-weighting estimator is invalid when treatment affects the outcome within a time interval, and we propose a modified version of the method that restores its validity in this setting. Finally, we analyze the consequences of choosing time intervals that are either excessively wide or excessively narrow, thereby formally establishing the need for time partitioning and providing practical guidance for selecting an appropriate partition based on the clinical setting.

What determines the unavoidable sample cost of learning cyclic causal structure? For cyclic linear non-Gaussian models, we study exact condensation recovery from observational data: identifying the strongly connected component (SCC) partition and all edges between components. We establish the first information-theoretic lower bounds on sample complexity for this target. For $p$ variables, maximum SCC size $s_{\max}$, and maximum external-parent count $d_B$, any estimator requires order $s_{\max}\log(ep/s_{\max})+d_B\log(ep/d_B)$ samples in the worst case over a regular model class. These bounds distinguish the costs of SCC membership and external-parent selection. Under principal invertibility and without correlation faithfulness, we establish a population block-exogeneity principle that identifies unknown root SCCs through residual independence and inclusion minimality. A sparse-adjustment characterization shows that small adjustment sets suffice to identify SCCs and their direct external parents, without regressing on all previously recovered variables. These characterizations yield BlockExo, which attains a structurally matching sample bound without knowing $s_{\max}$ or $d_B$ under suitable conditions. Simulations support the structural dependence of our sample bound and demonstrate BlockExo's sample-efficient recovery in comparisons with other methods for cyclic causal discovery.

We introduce the Functional Spectral-Newton Method (FSNM) for learning the leading singular structure of a conditional expectation operator without fixing a basis or reproducing kernel Hilbert space. FSNM fits a low-rank representation of the centered joint-to-product density ratio kernel by alternating functional Newton updates. Each update reduces to a preconditioned regression, which we approximate with vector-valued regression trees in a stagewise boosting procedure. At the population level, we establish descent and an $O(1/T)$ best-iterate block-stationarity rate under a relative weak-learner accuracy condition, and show that every nondegenerate local minimum over the full centered $L^2$ spaces is a globally optimal rank-$d$ approximation. Synthetic experiments show that FSNM recovers a low-rank density ratio and its leading spectral structure, and that the same learned kernel can answer multiple conditional queries without refitting.

Motivated by blind denoising in diffusion models, we study estimation of an unknown Gaussian noise level from a single high-dimensional observation, assuming the signal law P is known. We characterize the minimax mean-squared error under two structural assumptions on P. For signals with covering complexity k, the minimax rate is $\widetilde{\Theta}_\Lambda(\min\{\Delta_\Lambda^2,d^{-1}+k^2d^{-2}\})$, and the constrained MLE attains it up to logarithmic factors. For $\alpha$-strongly log-concave signals, the rate is $\Theta_\Lambda(\min\{\Delta_\Lambda^2,(1+\alpha^{-1})^2d^{-1}\})$, attained up to constants by a $P$-centered norm estimator. These results show that the structure of the signal law determines both the difficulty of blind noise estimation and the appropriate estimator.

Multiobjective Bayesian optimization (MOBO) with Gaussian process (GP) surrogates is a sample efficient approach to solving multiobjective optimization problems. In MOBO, a Bayesian decision theoretic acquisition function guides the adaptive selection of new candidate inputs, on which objectives and constraints are evaluated to update the surrogate model sequentially. Existing approaches maintain independent GP models for the objectives and constraints, with new observations evaluating all objectives and constraints in a coupled fashion. However, the objectives and constraints often contain inherent correlations which, if exploited, can enable ${decoupled}$ evaluations where only a subset of them are evaluated at each round. We present a new approach that leverages a multitask GP model to jointly learn all objectives and constraints, and propose a total correlation metric that enables identifying an ${optimal}$ subset of objectives and constraints to be evaluated at every round, even under uniform evaluation costs. Theoretically, we show that our acquisition policy is asymptotically consistent despite decoupling and that our proposed decoupled subset selection rule maximizes the expected posterior entropy reduction about unevaluated tasks under mild conditions. Empirically, we show that our approach outperforms coupled and decoupled baselines in the state of the art.

Semi-competing risks occur when individuals may experience a non-terminal event and a terminal event, where the terminal event censors the non-terminal event but not vice versa. In the presence of covariate-dependent censoring, augmented inverse probability of censoring weighting (AIPCW) framework has not been developed for the special setting of semi-competing risks which, unlike competing risks, have asymmetric event time structure and complex estimands. Using semiparametric theory on coarsened data, we carefully develop an AIPCW framework for semi-competing risks. When treatment effects are of interest, we further integrate this approach with augmented inverse probability of treatment weighting (AIPTW), yielding a framework for estimating causal estimands under covariate-dependent censoring. The resulting estimators are shown to be doubly robust. Using the proposed framework, we estimate treatment specific risks of: i) non-terminal event, ii) terminal event without the non-terminal event, and iii) terminal event following the non-terminal event. We evaluate the finite sample performance through simulations, and apply the method to data from the Honolulu Asia Aging Study to assess the causal effects of midlife heavy drinking on late life cognitive impairment and mortality.

Goodness-of-fit tests for logistic regression are least reliable where they are most needed: at small samples their levels drift from the nominal one, and combining them worsens the drift. We propose a test whose statistic is a convolutional network, trained once on simulated departures, that reads misfit as a picture: a grid of standardized residuals over covariate ranks. The analyst never trains. The network ships frozen, and the p-value is the rank of the observed score within the analyst's own bootstrap, so the level is a property of the calibration rather than of what the network learned. We prove exactness under pivotality, asymptotic exactness without it, and a consistency theorem whose key condition is computable from the frozen weights in one forward pass, giving a per-alternative certificate; we also measure the test's blind cone. On a pre-declared sixty-cell grid the deployed level stays in the nominal band in fifty-eight cells, while that grid's power criterion failed; on a published benchmark it is the most stable of thirteen levels across settings and sample sizes, and the test outpowers every partition test at every sample size while ranking seventh overall. A bank-failure application and a versioned release close the paper.

We study distributionally robust reinforcement learning (DR-RL) in the average-reward setting under weak communication. Our main result provides finite-sample guarantees for estimating the robust optimal average reward and learning a near-optimal policy, covering both SA-rectangular and S-rectangular structures with divergence-based and distance-based uncertainty sets. Specifically, for Kullback--Leibler and $f_k$-divergence balls, we establish explicit radius conditions under which the robust average-reward Bellman equation admits a constant-gain solution, while for total variation and Wasserstein balls, any positive radius suffices without requiring the nominal MDP to be weakly communicating. Our algorithm is prior-knowledge-free and achieves sample complexities of $\widetilde O(|\mathcal{S}||\mathcal{A}|p_{\wedge}^{-1}\operatorname{Span}^{2}(u_{\delta}^{\ast})\epsilon^{-2})$ for estimating the robust optimal average reward and $\widetilde O(|\mathcal{S}||\mathcal{A}|p_{\wedge}^{-2}\operatorname{Span}^{2}(u_{\delta}^{\ast})\epsilon^{-2})$ for learning an $\epsilon$-optimal policy. Here, $p_{\wedge}$ is the smallest positive nominal transition probability and $u_{\delta}^{\ast}$ is a robust optimal bias function. We further provide an almost-tight explicit upper bound on $\operatorname{Span}(u_{\delta}^{\ast})$. Finally, we validate the predicted $n^{-1/2}$ convergence rate through numerical experiments.

From a deck of cards, how many cards do I need to draw, until the first ace? I identify the distribution for this waiting time $T$, and its satisfyingly nice expected value ${\rm E}\,T=(N+1)/(n+1)$, with $N$ the number of cards and $n$ the number of aces; hence $53/5=10.6$ for the standard setup. After having solved this Question One I go on to certain alternative solutions and extensions, involving e.g. Beta approximations. I also consider the distributions and means for the 2nd, the 3rd, the 4th occurrences of aces, with generalisations, where there is a Dirichlet distribution in wait for us, with further links to order statistics for the uniform. Furthermore, an apparatus is developed for obtaining estimators and full confidence distributions for applications where one knows the number $n$ of aces, but not the deck size $N$; and correspondingly for inference about the unknown population size $N$ when $n$ is known. If you have 1000 people in a room, and need to interview 11 of them until you've found the first left-handed person, how may left-handed are there in the room -- here we need both an estimate and a clear measure of uncertainty.

World models for epidemiology learn from records shaped by study designs, including assignment, sampling, measurement, and related processes. A model may therefore reconstruct observed trajectories while learning an intervention contrast that depends on how records were collected. We formalize study design as constraints on a world model latent-world interface, action mechanism, observation likelihood, and target readout. This distinguishes design-respecting models, which encode these constraints, from design-ignoring models, which fit selected records without relating assignment and observation to the intended intervention question. Using a large scale cluster-randomized test-negative trial, we hold the latent structure and fitting settings fixed and compare a design-ignoring case-count model with a design-respecting test-negative observation model. Both models achieve comparable factual reconstruction. Yet across 500 paired resampling experiments at different relative sampling intensity, median contrast changes are substantial for the design-ignoring model, versus merely marginal for the design-respecting model. The same qualitative separation holds when sampling also varies across clusters. Thus, factual reconstruction alone does not establish design alignment; testing whether fitted contrasts preserve design-implied observation-process invariances provides a sharper evaluation.

We study inference for hyperuniformity from one large realization of a stationary point process. For point processes with an essentially free translation action of completely positive entropy, empirical Fourier statistics at Lebesgue-almost every frequency have asymptotically complex Gaussian laws, without quantitative mixing or cumulant-summability assumptions. Randomly sampled frequencies then give a tractable limiting experiment for low-frequency spectral mass. If the frequencies are uniform on the Euclidean ball $B_r$ of radius $r$, the mean of the limiting squared Fourier statistic is \[ A_r=\frac{\sigma_\eta(B_r)}{\rho_\eta\lambda_d(B_r)}. \] Hyperuniformity is characterized by $A_r\to0$ as $r\downarrow0$, so no continuity assumption on the structure factor at the origin is needed. Under \[ A_r=s+c r^\alpha+O(r^\beta),\qquad 0<\alpha<\beta, \] with a specified bound on the remainder and local square-integrability of the normalized Bartlett density, a two-radius extrapolation eliminates the $r^\alpha$ term and estimates $s$ with an explicit $O(r^\beta)$ bias bound. This yields confidence bounds and a one-sided test of $H_0:s=0$ whose asymptotic rejection probability under each fixed null process does not exceed the chosen significance level.

Task performance need not determine which intervention mechanism an agent retains. We study causal retention: whether a frozen learned state answers a mechanism-probe map fixed independently of training, including action, context, direct target, value, and delay. For finite structural causal model classes, the optimal probe error is a Bayes decision risk. It vanishes exactly when every learning-interface fiber lies within one probe-answer fiber; any state obtained by post-processing that interface inherits the same lower bound. A posterior-coverage theorem characterizes budgeted retesting, while an exact edit decomposition shows that the shifted set is the unique support of an error-free target update. Causal Core implements these conditions through evidence-gated writing, readout filtering, temporal credit, hidden-context setup, and local diagnostic updates. Experiments cover finite causal systems, continuous simulators, an official TD-MPC2 world model, and Qwen2.5-7B-Instruct. A frozen Qwen last-layer probe reaches 0.958 balanced accuracy on source mechanisms but 0.583 on changed delays; the gated mechanism state reaches 1.000 and accepts only 0.056 of synchronized-readout candidates. In TD-MPC2, five target states per actuator recover effect-sign accuracy from 0.057 to 0.948 without degrading stable responses. Causal retention is therefore distinct from task sufficiency and source-domain decodability.

The purpose of this article is to explicitly describe how the integration technique known as the change of variables is used in graduate-level Probability and Statistics courses. Throughout the article, we employ this integration technique seven times across various examples and results. We demonstrate that four applications of this technique are required to derive the probability density function (PDF) of the Student's t-distribution.

The multivariate expectile-based distribution is a `spiked' perturbation of a multivariate Gaussian distribution, proposed by Arbel et al. in 2023. The rejection sampler proposed for this distribution in the initial work is valid, but its acceptance probability can degenerate badly both in high dimension and in the strongly asymmetric limit. We give a simple alternative. After affine whitening and a polar decomposition, the sampling problem reduces exactly to a univariate distribution, and an additional hyperbolic change of variables exhibits this distribution as log-concave, so that a universal and uniformly efficient construction of Devroye applies. The resulting exact sampler therefore enjoys an acceptance probability of at least $1 - \exp \left( -1 \right) = 0.632120 \ldots$, uniformly over the dimension and all admissible asymmetry parameters.

In observational causal inference, covariate balancing is widely used to reduce source-target covariate shift, but under weak overlap in high dimensions, stronger balance can induce concentrated weights and increase variance. Balance measures how well the target covariate distribution is represented, but does not by itself determine how reliably the counterfactual mean can be estimated. We develop risk-calibrated balancing for the average treatment effect on the treated, which applies ridge augmentation to any normalised base weights and selects its penalty using conditional prediction risk of the counterfactual mean. Under a random-effects predictive model, we derive an exact finite-sample decomposition of this risk into residual covariate imbalance and weight-induced variance. For design-independent base weights under proportional asymptotics, we characterise how limiting risk depends on source and target covariance geometry, population mean shift, and weight concentration. For covariate-adaptive base weights, we develop a uniformly consistent target-aware risk estimator whose minimiser attains vanishing scaled oracle excess risk. Simulations show that the high-dimensional risk predictions remain informative for adaptive balancing and that target-aware tuning generally reduces excess target risk. Empirical analyses of job-training and single-cell perturbation data show that risk-calibrated balancing generally improves on the corresponding base estimators, with larger gains under weaker overlap.

We introduce a network realized GARCH-It\^o model in which volatility transmission is a dynamic relation among the latent daily integrated volatilities of multiple assets. An unknown directed and signed network is embedded in a continuous-time variance process and appears in the resulting exponential daily recursion. Intraday returns identify integrated volatility and hence the network that governs its propagation. For a fixed network template, feasible estimation based on realized volatility is first-order equivalent to estimation based on latent integrated volatility. For an unknown network, regularization selects the relevant structure, and, under oracle conditions, a fixed-rank local refit provides conditional inference for model-implied response-based connectedness. In an application to nine U.S. sector ETFs from 2007 to 2025, the model attains the lowest average out-of-sample QLIKE among the reported recursive forecasts, although its advantage over HAR is not statistically significant. In the full sample the weighted LASSO selects an empty sparse support, so inference concerns the rank-one projected factor component; this conditional analysis identifies Energy as a net volatility transmitter. Aggregate connectedness is more stable across volatility measures than individual sparse channels.

We study approximate sampling: given $N$ independent samples from a proposal distribution $\mu$, the goal is to select one whose distribution is close to a target $\pi$ specified only up to a normalizing constant. Block and Polyanskiy (2023) provide finite budget error bounds for approximate rejection sampling (RS) as a function of the acceptance threshold $M$. The threshold $M$ giving the smallest bound, however, depends on properties of $(\pi,\mu)$ that are typically unavailable from the observed sample. This raises a natural question: Can one attain the best RS guarantee without taking $M$ as input? We answer affirmatively by proposing a parameter-free sampling algorithm called uniform race (UR), based on importance weights, which are ratios of target to proposal probabilities (or densities). It divides each observed weight by an independent uniform random variable to form a score and returns the candidate with the largest score. For every budget $N$, its total variation error satisfies the RS upper bound for every fixed threshold $M$ simultaneously, thereby achieving the best such bound in hindsight. We also characterize its output distribution conditional on the largest score, identifying when it is exactly the target $\pi$. Uniform race has no larger total variation error than a natural budget-calibrated RS derived from Rohatgi et al. (2025) and sampling importance resampling (SIR). In particular, we exhibit instances where UR's error is exponentially smaller in $N$ than that of either baseline. Furthermore, we establish conditions under which attaining this RS guarantee for every $(\pi,\mu)$ uniquely determines the selection probabilities as those of UR. Finally, test-time scaling experiments on LLM math-reasoning tasks corroborate the theoretical comparisons and demonstrate that UR remains competitive in ground-truth accuracy without requiring threshold selection.

Statistical model transformations based on randomly stopped sums, maxima and minima are widely used to extend statistical models. We characterize the complete set of stopping models for which randomly stopped sum and extreme model transformations function as statistically stable (idempotent) model extensions. Stability requires the underlying stopping model to be closed under pgf composition. We prove that any finite-dimensional, connected stopping model closed under pgf composition is necessarily a family of random variables whose pgfs commute. Using the corresponding Koenigs function, we establish that these models form a statistical manifold admitting a global, one-dimensional parametrization $\theta = \Pr(N=1) \in (0, \theta_*]$, where the probability mass at $i$ is a polynomial in $\theta$ of degree at most $i$. Finally, we establish a duality between stopping models closed and containing the identity variable (the ones yielding stable extensions) and the set of probability distributions supported on the positive integers. These findings disprove the long standing conjecture that statistical stability occurs only under geometric stopping.

Estimating conditional statistics and learning representations of a population of conditional distributions are central problems in many data-driven applications, including uncertainty quantification and dynamical systems analysis. Conditional mean operators (CMOs), a class of linear operators between function spaces, resolve these objectives by providing access to a broad class of conditional statistics. However, existing methods typically estimate each CMO independently or constrain it to prespecified function spaces, thereby preventing the exploitation of shared structure across related distributions. In this work, we posit that related CMOs share finite-dimensional input and output function spaces, and are specialized for each task with a linear operator mapping these spaces. Based on this hypothesis, we introduce MTL-CMO, a multi-task framework that jointly learns shared function spaces and task-specific operators across multiple datasets. We further introduce T-CMO, a transfer learning method that reuses the shared spaces to estimate, in closed form, the operator of a new conditional distribution. We establish statistical guarantees quantifying the benefits of jointly learning the shared function spaces. Our experiments demonstrate that learning shared function spaces improves uncertainty quantification across a broad range of conditional distributions and, when applied to Langevin and plasma dynamics, yields compact representations of complex dynamics that retain physically meaningful information and enable parameter identification.

We establish a Donsker-Prokhorov invariance principle in some Besov spaces. Specifically, we show that polygonal line processes associated with partial sum processes of a triangular array of row-wise independent random variables converge in distribution to Brownian motion, extending earlier results in the literature.

Nested sampling estimates the evidence by assigning prior volumes to an ordered sequence of likelihood contours. Markov chain constrained-prior samplers, which are used in high-dimensional problems, generate many intermediate states before producing the next classic sample. These intermediate states are discarded because their correlation prevents them from being inserted into the ordered NS sequence without changing its order-statistic law. This paper introduces a novel method of using them to improve evidence estimation, by formulating NS in a Bayesian way and conditioning on phantom samples as Monte Carlo observations. We then present the open-source software package, JAXNS v3, and its implementation choices. We validate the approach on a set of problems, and identify its limitations via ablation. In our experiments, when problem structure is well resolved, conditioning on all retained phantom samples reduces log-evidence RMSE by at least $30\%$, with larger improvements at higher dimensionality. For the tested problems with unresolved structure, full phantom conditioning produces no detectable improvement or deterioration in evidence accuracy. Phantom conditioning produces overconfident evidence uncertainties. We also introduce two dynamic nested sampling allocation schemes. Evidence-improving allocation approximately halves the number of likelihood evaluations required to achieve comparable evidence accuracy relative to uniform allocation. Posterior-improving allocation doubles the classic posterior's Kish effective sample size for $17.5\%$ additional likelihood evaluations.

We here develop a functional neural network, termed MatFAE, for learning trajectories on the Riemannian manifold of symmetric positive definite (SPD) matrices. MatFAE features intrinsic layers that map manifold-valued functions to Euclidean vector-valued functions, followed by a functional layer that projects them into a finite-dimensional Euclidean space. Unlike most neural networks for discrete-time sequences, MatFAE treats each sequence as a continuous function and can therefore encode trajectory dynamics (e.g., first-order derivatives) in its latent representations. Additionally, the morphology of the functional weights in the functional layer offers interpretability by revealing the regions of the input functional data that contribute most to the latent representations. We justify the design principles and properties of each intrinsic layer and detail how matrix factorization is handled during backpropagation. We apply MatFAE to a range of fMRI datasets, demonstrating its ability to efficiently learn informative representations from high-dimensional SPD trajectories and its practical value for real-world neuroimaging analysis.

Score-based diffusion models are increasingly considered in settings where the underlying data distribution may differ from the training distribution, yet existing theoretical guarantees largely focus on the no-shift setting. In this work, we study robust score estimation under Wasserstein perturbations of a reference distribution. For the Ornstein--Uhlenbeck diffusion, we show that robust estimation decomposes into two fundamental components: the statistical cost of learning the reference distribution and the intrinsic cost of distribution shift. The latter scales quadratically with the Wasserstein radius, and this dependence is minimax optimal. We construct an explicit finite-sample estimator achieving the resulting robust minimax rate without knowing the shift radius. When the reference distribution lies on an unknown low-dimensional subspace, the statistical term adapts to the intrinsic dimension while the shift cost remains unchanged. Finally, we show that the same decomposition governs positive-time reverse sampling and obtain matching minimax guarantees in KL divergence. Together, these results characterize how finite data, intrinsic dimension, and distribution shift affect the robustness of score-based diffusion models.

The standard approach to stochastic simulation from the von Mises-Fisher distribution is a rejection sampler proposed by Ulrich and Wood. This note provides a theoretical justification of its efficiency, lower-bounding its acceptance probability uniformly over the dimension and concentration parameter. A novel interpretation of the proposal is given which provides some intuition for this efficiency.

We study the nonparametric statistical inverse problem of recovering the space-dependent permittivity in a coupled Fokker-Planck-Darcy system from discrete, noisy observations of the Fokker-Planck solution. We consider the coupled parabolic-elliptic system on a bounded domain with the physically natural no-flux boundary condition for the Fokker-Planck equation and homogeneous Dirichlet boundary conditions for the Darcy equation. The inverse problem is indirect, since the unknown coefficient enters only through the elliptic equation and is observed only through its effect on the density in the parabolic equation. We first develop the analytical theory of the forward problem required for the statistical analysis, establishing well-posedness and a priori bounds uniform over the admissible set of permittivities. We then derive two stability estimates for the inverse problem, namely a Lipschitz-type forward estimate and a generalized backward estimate. Placing a rescaled Gaussian process prior on the log-shifted permittivity, we show that the posterior contracts around the truth at an explicit polynomial rate in the number of observations, with the posterior mean converging at the same rate. Numerical experiments in one and two dimensions, using preconditioned Crank-Nicolson and ensemble Kalman filter algorithms, complement our theoretical results.

Groups differing in social and economic circumstances often differ markedly in depressive symptom levels, yet their estimated symptom networks show few clear differences. Such null results are often read as evidence that the groups' symptom systems are the same, but they only show that no difference was detected in how symptoms interact. How readily symptoms become present, and how this depends on the contextual conditions under which they are observed, is rarely compared. Building on a slow-fast perspective, we treat depressive symptoms as a relatively fast-changing system embedded in more persistent social, economic, psychosocial, health, and lifestyle conditions that shape symptom activation. Using a cohort of 23,689 adults from Amsterdam, we illustrate this with PHQ-9 depressive symptoms and a composite index of persistent socioeconomic, psychosocial, and health-related conditions, which we call Slow Risk Load (SRL). Participants with high SRL had substantially higher symptom levels than those with low SRL. A standard Network Comparison Test detected no overall difference in network structure, although global strength was modestly higher in the high-SRL group. Ising models, which separate how symptoms co-occur from how readily each becomes present, indicate higher activation under high SRL for eight of nine symptoms, and model comparison consistently favored group-specific activation over group-specific interactions. The slow-fast perspective reframes group network comparisons as a multiple-timescale problem in which network structure, symptom activation, and context should be considered jointly. It places symptom-focused and context-focused interventions at different layers of the same coupled system rather than in competition.

AI-workflow governance cannot be reduced to checking the final answer: an apparently safe answer may rest on a fragile evidence path that ordinary evaluation cannot see, localize or govern. We call this hidden fragility a "structural iceberg": hallucinations and unsupported claims may form its visible tip, while consequential weakness remains submerged. Stochastic semantic evidence graphs (SSEGs) expose these icebergs by preserving workflow channels, propagating local uncertainty and identifying the hidden paths on which an apparently safe output depends. ALCE and RAGTruth show that visible failures at the tip -unsupported citations and hallucinated spans -rest on distinct submerged weaknesses and therefore require different interventions. Across retrieval, tool-use and controlled stress tests, SSEG localizes those weaknesses, supports targeted repair, produces no false automatic passes in 35,000 known-truth cases and reduces ToolSandbox review by 28.8% across 96 executions from two agent models. The same structural view carries into end-to-end governance: in a separately sealed 1,200-case FinGovBench study, adding SSEG to GPT-OSS-20B reduces unsafe releases from 452/660 to 8/660 while releasing all 540 safe cases and correctly distinguishing 592/600 matched workflow pairs. An unchanged-gate transfer to Qwen3-8B releases all 540 safe cases and none of 660 unsafe cases, whereas flat-UQ releases 520 unsafe cases. SSEG therefore moves governance below surface-level output checking, turning hidden evidence dependencies into auditable, path-specific decisions about intervention, revalidation and release.

Transformers have emerged as powerful architectures for learning solution operators of physical systems. Empirically the prediction error has been observed to decrease when the data size and model size increase, suggesting neural scaling behavior. Yet a theoretical understanding of such scaling laws for transformer-based operator learning remains limited. In this work, we develop a theoretical framework for characterizing the approximation and generalization errors of transformer-based operator learning. Our analysis builds on a local-to-global approximation principle that is naturally aligned with the softmax attention mechanism and yields discretization-invariant output functions. On approximation theory, we derive a universal approximation error of transformer-based operator learning for H\"older-regular operators. On generalization theory, we establish a power scaling law between the prediction error and the training data size. The rate of convergence represented by the scaling exponent explicitly reflects the dimensions of the input and output domains, the regularity of the underlying functions and operators, and crucially, the intrinsic dimension of the input function class. By exploiting this intrinsic low-dimensional structure, our analysis yields a power-law generalization rate for operator learning, in contrast to the logarithmic-type power-law rates appearing in existing analyses of operator learning with feedforward neural networks. Numerical experiments validate the predicted power-law scaling and confirm that the convergence rate varies systematically with the intrinsic dimension of the input function class.

We introduce RareTrap, a framework for estimating the probability of severe behaviors in black box large language models (LLMs). A key challenge for probability estimation is defining a tractable distribution over the input space. To accomplish that, RareTrap uses a surrogate LLM and constructs a geometry-aware mapping from a lower-dimensional latent reference space into its token-embedding space to induce an explicit and reproducible distribution over input prompts. A response-level performance function is utilized on the response to quantify behavior severity. This enables sequential rare event simulation that concentrates evaluations on progressively more severe behaviors while preserving probability under the induced prompt distribution, which would otherwise be prohibitive to measure. Across 10 open-weight and two frontier models (GPT-5.4 and Claude Sonnet 4.6), we find that RareTrap successfully induces severe resource consumption behaviors and computes their probability with as few as 200 evaluations. RareTrap provides model developers a principled approach for evaluating language models under a common distribution, and prioritizing alignment effort to improve safety and mitigate risks.

Generating independent, equilibrium samples of molecular systems at scale remains a central obstacle in computational statistical mechanics. Boltzmann Generators address this by pairing a generative model with importance sampling to obtain consistent samples from the target distribution. We introduce Normalizing Flow Flow Maps (NF$^2$M), which combines the strengths of recent stochastic flow maps with the tractability of classic normalizing flows to build a Boltzmann Generator. Unlike most methods, which correct the generative model only at the end, NF$^2$M reweighs each denoising transition as generation proceeds, avoiding wasted compute on trajectories that are ultimately discarded. At each denoising step, a conditional normalizing flow proposes clean configurations given the current noisy state (a simpler task than sampling directly from the target) and its exact likelihood enables correcting each proposal toward the true denoising transition of the target Boltzmann distribution. This is in contrast to most existing methods, whose likelihoods are approximate or expensive to evaluate, undermining the statistical reliability of the correction. We establish consistency of the corrected transitions and bound how approximation errors propagate through the sampling chain. We evaluate NF$^2$M on peptide systems, demonstrating improved sampling efficiency and sample quality.

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

We study the high-dimensional recovery of a signal vector $\mathbf{x}$ in the presence of sparse Wishart-like noise. We define an $N \times N$ matrix $A = J+(\theta/N)\mathbf{xx}^{\top}$, where $\mathbf{xx}^{\top}$ is the rank-one deformation of the random noise matrix $J$. We consider a Wishart-like matrix $J={X}^{\top} X$, where $X$ is a sparse $M \times N$ random matrix with entries $X_{ij} = c_{ij}W_{ij}$, with $c_{ij}$ regulating the density of non-zero elements, and $W_{ij}$ the bond weights. Using the replica method, we compute analytically the top eigenpair statistics of $A$, and their dependence on the signal strength $\theta$, the rectangularity ratio $\alpha=\sqrt{M/N}$, and the average connectivity of the noise. The spectral observables are expressed in terms of a system of Recursive Distributional Equations, which are efficiently solved via a Population Dynamics algorithm. They allow us to compute the average largest eigenvalue $\langle\lambda_1\rangle_{A}$, the average top eigenvector component density, and the average overlap between the top eigenvector of $A$ and $\mathbf{x}$. We identify a critical threshold $\theta_{\mathrm{crit}}$--depending on the average connectivity of the noise--that marks a BBP-like phase transition: below this value, $\langle\lambda_1\rangle_{A}$ is unaffected by the signal, and the overlap vanishes. Thus, the signal is not recoverable from the top eigenvector of $A$. For $\theta>\theta_{\mathrm{crit}}$, the signal-related outlier eigenvalue becomes $\langle\lambda_1\rangle_{A}$ and the overlap is nonzero, allowing for recovery of the signal. The results are in excellent agreement with numerical diagonalisation. We show that in the dense limit, the recovery threshold and eigen-statistics converge to the results predicted by the classical BBP transition for additive rank-one deformations of dense Wishart matrices.

In ecology, psychometrics, and the analysis of social and financial networks, binary matrices are often analyzed conditional on their observed row and column sums, which restricts the problem to a finite sample space of matrices with the same margins. Two fundamental problems are to count this space and to sample uniformly from it. Sequential importance sampling (SIS) addresses both with independent weighted samples and an unbiased count estimator, but its efficiency depends critically on the proposal distribution. Existing proposals are analytically designed, and their accuracy can vary substantially with the margins. We show that the ideal SIS proposal, under which every weight equals the count and the variance vanishes, is exactly the policy of a generative flow network (GFlowNet) with unit reward on every matrix that has the given margins. We therefore propose MarginFlow, a framework that turns the design of the proposal into a learning problem and amortizes it across margins by exploiting their self-similarity. Every partial matrix is itself an instance with reduced margins, so one set transformer that reads the remaining margins serves every margin. We train MarginFlow on a pool of 1904 margins and evaluate it zero-shot on 1190 held-out margins, synthetic and real, from $3\times3$ to $870\times6$. On 1187 of the 1190 margins it matches or beats the best of 31 analytically designed configurations, chosen post hoc for each margin, and its median effective sample fraction is 99.8%. On the 56 margins where that best loses more than one nat of effective sample size, MarginFlow wins every one and raises the median effective sample fraction from 10.3% to 94.1%.

Plug-and-Play (PnP) methods replace proximal operators with learned denoisers, which produce state-of-the-art reconstruction quality, but sacrifice the variational interpretation and convergence guarantees of proximal methods. Learned Proximal Networks (LPNs) address this problem by designing the denoiser architecture so that it is exactly the proximal operator of a regularizer. In this work, we clarify the theoretical foundations of LPNs and extend the framework to a broader class of activation functions. We then study two practical mechanisms for controlling the learned prior: (i) averaging the LPN with the identity, a common heuristic in PnP methods, and (ii) directly scaling the implicit regularizer induced by the proximal operator. In particular, we characterize the regularizer induced by averaging and develop a convergent method for evaluating the scaled proximal operator.

Low-rank matrix and tensor bandits exploit structured interactions but typically assume a known reward link. Recent single-index bandit methods accommodate unknown links without directly exploiting matrix or tensor rank. We address this gap by studying stochastic matrix and tensor bandits with an unknown shared Lipschitz link and a low-rank index parameter under known regular candidate distributions and finite-variance noise. For monotone links, T-ESTOR combines robust, rank-adaptive Stein estimation with epoch-based greedy selection. Under exact selected-score access and a uniformly positive selected-design Stein signal, it achieves square-root regret with dimension dependence determined by the low-rank structure. For every admissible design, the monotone lower bound matches the rank, dimension, and horizon dependence up to logarithmic factors at large horizons, for fixed menu size and model/design constants. For nonmonotone links under a nonzero base-law Stein signal, T-BSTOR combines structured estimation with robust bin-based learning and attains the optimal $\widetilde{O}(T^{2/3})$ horizon rate for fixed dimensions, menu size, and model/design constants. Synthetic and CCLE-based experiments illustrate the benefits of structured estimation relative to vectorized and competing single-index baseline methods.

Estimating causal effects is challenging with observational data due to the lack of exchangeability, and is further complicated with spatial data. In real-life applications such as agricultural and environmental studies, unobserved spatial factors (spatial confounding, SC) and interactions of units in nearby locations (spatial interference, SI) commonly occur jointly, yet existing parametric methods address them separately. We develop a unified framework that treats both simultaneously. We first define the estimand precisely and give an identification result. We then study two estimators. The matching estimator iDAPS integrates neighbourhood exposure and spatial proximity into propensity score matching with balance-optimised weights. We prove a variogram bias bound showing its confounding bias is controlled by the spatial variability of the confounder at the matched distance; this justifies the composite metric, yields a consistency result, and provides a computable diagnostic. The doubly robust estimator recoverU+ augments the propensity and outcome models with the recovered confounder and the exposure; we prove it is doubly robust up to an explicit residual bias from the unrecoverable confounder component. We show that independence-based standard errors are anticonservative under spatial dependence; we provide spatial heteroskedasticity and autocorrelation consistent (HAC), block-bootstrap, and randomisation-based alternatives whose calibration we verify. Every theoretical result is validated numerically. Estimating the effect of selective catalytic and non-catalytic reduction (SCR/SNCR) technologies on ambient ozone reveals no evidence of ozone reduction, while a naive analysis would have delivered a confidently wrong sign. All methods are implemented in the developed open-source R package spaci.

Prediction models can perform poorly when the deployment population differs from the training population. Data from the target population would help, but individual-level target data may be inaccessible because of access restrictions or reporting conventions. We consider a multi-resolution setting in which individual-level data are available from a source population, while the target population is observed only through subgroup summaries. We propose SAGE, a one-step estimator that updates a source-trained predictor using a gradient estimated from these summaries. Motivated by diagnostics consistent with the random distribution shift model, we choose SAGE's step-size to account for both sampling and distributional uncertainty. Under this model, SAGE reduces mean asymptotic target excess risk relative to the source-trained predictor. We also show that, under the model, entropy-balancing weighted empirical risk minimization (EB), which reweights source observations to match the target summaries, is asymptotically equivalent to a full-step SAGE update. SAGE with the optimal step-size has asymptotic mean squared error no larger than that of EB. Across real-world datasets, SAGE generally improves on the source-trained predictor and one-step updates that ignore distribution shift, including when the random shift model only partially captures the observed shifts. Compared to EB, SAGE improves prediction more consistently across the sample size and shift settings studied.

LLM leaderboard gains can reflect selection among privately evaluated model variants, yet neither the number of variants nor their dependence is public. We ask how many hidden variants a published margin can support while retaining statistical evidence of a provider's advantage over a fixed comparator. For a fixed candidate family under a Gaussian margin model, we derive a sensitivity curve that reports this maximum count as a function of a lower bound on within-family correlation. The relevant correlation must match the score used for ranking and the sampling model: in a controlled family, pooled item correlation is 0.90, whereas composite-score correlation is 0.46 under item resampling and 0.92 when MMLU subjects are resampled. An item-based audit of 394 adjacent-rank claims on the Open LLM Leaderboard finds that 391 lack statistical support even before accounting for selection. Among claims that pass the uncorrected test, certification can depend on assumptions about the hidden family's correlation. The resulting curves make these assumptions explicit without estimating the unobserved search size.

Generalized Langevin equations describe non-Markovian dynamics in which the evolution of resolved variables depends on their past. We propose a memory-conditioned diffusion method for learning stochastic flow maps of these dynamics from observed trajectories, without identifying a memory kernel or reconstructing unresolved variables. A compact, recursively updated bank of exponential filters enables the flow map to retain predictive history over multiple time scales without conditioning on long observation windows. The next-step distribution is conditioned on the current observation and this memory state, whose storage and update costs are independent of the history length for a fixed bank size. Predictive criteria guide the memory budget, with reference-assisted selection in the vector benchmark, and an optional linear projection further reduces the conditioning dimension. A kernel-based score estimator generates conditional samples without training a score network, and these samples are used to train a neural flow map for autoregressive simulation. Three numerical examples assess long-memory retention at small conditioning dimension, predictive compression in coupled vector dynamics, and non-Gaussian conditional distributions and intermittent events. The non-Gaussian example reproduces conditional asymmetry and burst statistics in a stochastic model of the plasma scrape-off layer.

Approximate message passing on factor graphs underlies two dominant families of probabilistic inference algorithms: expectation propagation (EP) and variational message passing (VMP). Both methods approximate the marginal at each factor edge, forcing an iterative round-robin schedule, risking negative-precision messages, and, for VMP, collapsing to point estimates at Dirac-delta factors. We introduce Direct Message Approximation (DMA), which approximates factor-to-variable messages directly rather than the marginal. For normalisable factors, we define a consistency condition (requiring exactness when all other incoming messages are Dirac deltas) to guide message construction. We prove a master theorem (proper messages, any graph) bounding marginal KL from message KL, with three structural corollaries: Dirac-input consistency, no EP-style inner-loop iteration, and no negative-precision messages. Further, we prove a complementary $O(1/r^2)$ guarantee for the inherently improper backward message of the product factor, whose closed-form treatment has resisted prior work. As a concrete instantiation, we derive explicit DMA messages for the product and leaky-ReLU factors and assemble a Bayesian neural network (BNN) inference algorithm with one forward/backward sweep per training example and no gradient learning-rate hyperparameter, validating that the structural guarantees translate to predictive uncertainty that widens in data-sparse regions, including under model mismatch.

Recovering probability measures from their moments has numerous applications, esp. in connection with the method of moments in statistics and optimization. In the setting where measure need not be finitely atomic, but its support is known to be compact and semialgebraic with codimension at least one, the problem is still open. We combine moment-matrix kernel information with the Christoffel--Darboux kernel to provide a discrete approximation of the support. To validate the proposed approach, we test our algorithm on analytically computed moments and pseudo-moments arising from polynomial optimization problems without unique global minimizers. This complements well-known recent work on recovery of measures with algebraic support, where the kernel of a moment matrix can reveal polynomials vanishing on the support, and on recovery of sufficiently regular full-dimensional supports, where estimators constructed by thresholding the Christoffel--Darboux kernel are known to converge asymptotically to the support.

Estimation of parameter of probabilistic models is an important task in the field of machine learning.For models of discrete variables, calculation of the normalization constant of model is sometimes difficult and a lot of researches have been done to avoid the calculation of the normalization constant. In this paper, we tackle with the difficulty by combining a technique of empirical localization and a deformed Bregman divergence.The technique of empirical localization makes it possible to drastically reduce computational cost of the calculation of the normalization constant, and in addition, appropriate choice of the deformation for the Bregman divergence can invest the proposed estimator with various kinds of favorable statistical properties, such as efficiency or robustness against outlier noise.

Diffusion models have revolutionized generative modeling for continuous data through the gradual refinement of a belief state. This iterative refinement has not yet carried over to discrete diffusion models, which discard uncertainty at intermediate steps through categorical sampling (information collapse). We propose Simplex Diffusion Models (SDMs), a framework that lifts the diffusion process to the probability simplex to represent beliefs over categories. SDMs admit probability paths with closed-form reverse transitions and can be trained with a simple cross-entropy loss. Contrary to earlier proposals such as Dirichlet Flow Matching which requires integrating an ordinary differential equation, we introduce a DDIM-like sampler with a tunable level of stochasticity. Because SDMs operate on samples on the simplex, they can carry uncertainty across denoising steps, which mitigates information collapse. On OpenWebText, SDMs are competitive with strong Discrete Diffusion baselines, achieving $17.0$ GenPPL at $5.46$ unigram entropy in 64 sampling steps, close to real validation data. Even without Self-Conditioning (SC), SDMs outperform masked and uniform diffusion (with SC or predictor-corrector sampling) on code generation (TinyGSM, $T=0.1$; $49.0\%$ vs. $45.8\%$). Distilled down to 8 steps, SDMs solve $32.1\%$ of GSM8K problems, more than distilled Discrete Diffusion models with 128 steps ($21.4\%$).

We study the generalization of spectral gradient descent (SpecGD) in overparameterized matrix classification with corrupted labels. Each input combines a shared low-rank signal with a rank-one sample-specific perturbation, referred to as a shortcut, that enables memorization but does not generalize. We contrast collapsed shortcuts, which share a singular direction, with dispersed shortcuts, which occupy distinct singular directions. Changing only this geometry can reverse the relative generalization of GD and SpecGD: collapsed shortcuts can favor SpecGD, while dispersed shortcuts can favor GD. In the dispersed regime, exact shortcut orthogonality eliminates the signal from the late-stage SpecGD direction, while vanishing random correlations collectively generate a small but generalization-relevant signal through a second-order effect. To identify the direction selected by SpecGD, which the spectral max-margin problem alone does not determine, we combine a refined analysis of its dual with the exponentiated-gradient dynamics of normalized loss weights. Finally, we show that a single SpecGD step can already interpolate and generalize well, while continued training converges to a direction with substantially worse generalization.

We study the problem of learning directed acyclic graphs (DAGs) with nonnegative edge weights from observational data. We propose the Nonnegative and Concomitant (NoCo) DAG estimator, which jointly recovers the weighted graph structure and the exogenous noise variances in the linear structural equation model for the observations. Different from prior art, this noise-adaptive formulation blends a smoothed concomitant lasso criterion with a simpler log-determinant acyclicity constraint that exploits nonnegativity and yields a more benign optimization landscape. Specifically, nonnegative weights allow us to impose acyclicity directly on the adjacency matrix without elementwise squaring of its entries, thus avoiding the well-documented degeneracy of the Karush-Kuhn-Tucker conditions. Computationally, we develop a method of multipliers' algorithm that accommodates both homoscedastic and heteroscedastic noise profiles. Within each iteration, we use block successive convex approximation to minimize the augmented Lagrangian, alternating between proximal gradient steps for the adjacency matrix and closed-form updates for the noise scales. Simulated experiments demonstrate NoCo's improved structural and edge-weight recovery relative to competing methods in a variety of settings, highlighting the benefits of exploiting nonnegativity along with noise adaptivity.

Bi-invariant geodesic regression generalizes linear regression to Lie groups. Its main feature is that it respects the symmetries of the group so that the resulting estimator is independent of arbitrary choices such as a reference frame. However, the local existence and uniqueness of the underlying estimator have not been shown until now. Furthermore, the convergence properties of the proposed algorithm for computing the estimator are not known. In this work, we investigate these questions. We prove that, locally, a unique estimator exists and give explicit bounds on the size of this neighborhood. We also show that the proposed iterative algorithm converges linearly to this estimator.

Random forests (RF) are tree-based models that capture complex, high-order predictor interactions. In biomedical studies, ordered outcomes often have outliers, skewness, or mixtures of continuous and ordinal values (e.g., due to detection limits). Standard RF splitting targets conditional means via squared errors and is sensitive to such irregularities. We propose a Wilcoxon regression tree that selects splits by maximizing a rank-based impurity reduction, equivalent to maximizing the squared Wilcoxon rank-sum statistic. Because the criterion depends only on outcome ranks, it is invariant to monotone transformations and naturally accommodates continuous, ordinal, or mixed outcomes. We develop the Wilcoxon random forest (WRF), which aggregates Wilcoxon regression trees via subsampling and estimates conditional distributions using forest-weighted empirical CDFs. We establish consistency of the WRF distribution estimator under regularity conditions. We evaluate distributional prediction using calibration diagnostics and continuous ranked probability scores, and define an out-of-bag permutation variable-importance measure. Simulations show that the WRF performs comparably to the standard quantile regression forest when errors are symmetric and homoscedastic, and yields improved quantile estimation and calibration when outcomes are skewed or heteroscedastic. We apply the WRF to predict CD4 cell count and HIV viral load six months after antiretroviral therapy initiation in a multicenter Latin American cohort, and the WRF improves threshold probability estimation and accommodates the large mass at the viral load detection limit.

A learner's current observations need not determine its response to further training. We formulate intrinsic-extrinsic coupling through the continuation-conditioned value of a constrained learning-state intervention, with observation-relative fibers describing present agreement. An executable finite-frame classifier-head write protects current logits while repairing specified historical margins under finite-precision acceptance checks. We distinguish local admissibility, continuation-conditioned intervention value, and complete-policy performance. A matched four-cell contrast identifies readout-specific non-additivity between the same intrinsic intervention and alternative external continuations. In a CLINC-derived class-incremental setting, replay changes the write's 32-update contribution from five correct predictions to zero. Nonzero interactions also occur under output distillation, with a RoBERTa backbone, and under optimizer-native SGDW dynamics. Under SGDW, correct-count interactions are negative in all three activated roots at 128 updates, showing that coupling need not imply positive synergy. The mathematical analysis distinguishes feasible local repairs and favorable terminal outputs from training-reachable repair regions. Separate coordination tests show that content controls match or exceed the development gain, while a five-root fresh-test comparison with Fiber present in every arm shows root-dependent rather than uniformly beneficial correct-count effects. On the secondary cross-entropy readout, guided allocation yields lower mean loss than standard replay in all five pairs. Together, these results make intrinsic-extrinsic coupling operational by connecting executable state geometry to continuation-conditioned value, matched interaction identification, and closed-loop coordination, while separating identified coupling from complete-policy performance.

Large language models can interpret natural lan- guage, yet robust decisions remain challenging. Jev-like models expose structured choices, but these interfaces do not directly provide numeri- cal values at a requested precision. We propose NUMERICJEV, a training-free numerical decod- ing algorithm that enables numerical output from any LLM with a Jev-like structured-choice in- terface. Surprisingly, on our arithmetic bench- mark, it outperforms direct selection from a can- didate list containing the correct answer by 2.93 percentage points (Figure 1). Our motivation comes from the observation that numerical range selection is itself a decision problem that Jev- like LLMs can address. NUMERICJEV recur- sively refines a range through a multiway deci- sion tree while retaining the original question in context, without parameter updates or hidden- state access. On a 100-value grid, a ten-way tree requires only two decision rounds. Range- normalized MAE is 1.84% versus 5.18% for di- rect choice. A separate three-date historical- index study yields 4.58% mean relative recall er- ror and 0% readout error when the value is sup- plied. Code is available at https://github. com/Bring-AI/jev-numeric.

Inspired by the problem of identifying the best model from a collection of large language models (LLMs) with heterogeneous querying costs, we formulate and analyse a variant of the multi-armed bandit (MAB) with (i) dueling feedback, where pairwise comparisons between model responses provide robust preference signals, and (ii) heterogeneous sampling costs, reflecting the differing costs of querying different LLMs. Assuming the existence of a Condorcet winner, a condition we empirically validate across multiple real-world datasets, we propose a Track-and-Stop style algorithm for best-arm identification with prescribed confidence. We prove that the algorithm almost surely achieves the asymptotically optimal cost as the error tends to zero. Finally, we extensively evaluate our approach on both synthetic and real-world instances, demonstrating consistent improvements over classical cost-unaware algorithms and their cost-aware extensions.

LLM-as-a-judge has become the de facto standard for scalable, subjective evaluation, yet current leaderboards compensate for systematic measurement bias by running ever more comparisons, an approach that is both statistically unsound and computationally wasteful. The root cause is an incomplete measurement model, treating LLM judges as neutral, interchangeable instruments ignores documented biases like position bias, verbosity bias, judge severity, and self-enhancement, that no volume of additional data can eliminate. We propose a unified latent variable framework that jointly models pairwise and ordinal data while explicitly correcting for these confounders, recovering reliable rankings from substantially fewer comparisons. Because fitting this model costs negligible compute relative to a single round of LLM inference, bias correction is not only more statistically rigorous but also a more sustainable approach to trustworthy evaluation.

Linear mixed-effect panel data models are widely used in longitudinal biomedical, environmental, and social studies, but are often sensitive to data contamination and numerous covariates. To address these challenges, we propose a robust variable selection approach based on a least squares approximation (LSA) of the density power divergence (DPD) objective function combined with the Adaptive LASSO penalty. The LSA converts the nonlinear DPD objective into a computationally efficient quadratic approximation while preserving the robustness of DPD estimation. Under suitable regularity conditions, the proposed DPD Adaptive LASSO-LSA estimator is shown to possess oracle properties, including selection consistency and asymptotic normality. Simulation studies demonstrate improved robustness, better sparsity recovery, and significantly enhanced computational efficiency than penalized likelihood methods, while maintaining performance comparable to existing DPD-based approaches. An application to the SWAN bone mineral density dataset illustrates the practical relevance of the proposed methodology for robust estimation and reliable identification of important covariates.

We study whether benchmarks reliably guide model selection as developers adapt to evaluation feedback across multiple criteria. We find that the worst-case test-set size needed to estimate the best score among $k$ adaptively chosen models, under any convex combination of the criteria, grows exponentially with the number of criteria, reaching the $\Theta(\sqrt{k})$ cost of answering $k$ adaptive statistical queries with only $O(\log k)$ criteria, at fixed accuracy and confidence. In attacks on multi-task large language model benchmarks with five to ten criteria, feedback restricted to nondominated task profiles produces large reused-to-held-out score gaps and frequent false winners. These results challenge a prominent explanation for prior observed reliable benchmark reuse---that developers mainly respond to convincing improvements over the current best---in rich-feedback settings, while leaving open how often ordinary model development encounters this vulnerability.

Short-term vacation rentals, as promoted by platforms such as Airbnb, Homeaway, Vrbo, etc., are a growing component of the travel industry. This paper provides a unique, large dataset on global market occupancy for the short-term rental market using data from the American company Wheelhouse. The dataset consists of data for $500$ markets around the world. For each market, a daily occupancy time series from January $2017$ to December $2022$ is provided, allowing for studies of local and global patterns in the evolution of the short-term rental market. Additionally, the dataset includes curves representing the booking trajectory of each market and stay date up to one year prior to the stay date. This large dataset comprises a unique combination of time series and survival analysis data, and is suitable as a methodological benchmark for both classical statistical and machine learning models.