Papers

Evaluating the average causal effects of treatment strategies on recurrent event outcomes, such as heart attacks or renal failure, is important in clinical and medical research. However, the analysis becomes increasingly complex as multiple interacting factors are considered within a longitudinal setting. In this paper, we use advanced methodologies to estimate the average causal effects of standard versus intensive blood pressure-lowering therapies on acute kidney injury recurrences. We address time-varying treatment and confounding, and model misspecification during the identification and estimation processes for the effect estimands. We analyze the Systolic Blood Pressure Intervention Trial data set using our proposed method, accounting for medication adherence and the semi-competing risk of death observed in the data.

An Edgeworth expansion of first order is established for general linear rank statistics under the null hypothesis with a remainder term that is usually of order $n^{-1}$. Furthermore, corresponding results for the second order are formulated, but not proved here. The proof for the first order is based on Stein's method and on an extension of the combinatorial method of Bolthausen. It is also shown that conditions of van Zwet imply up to a small factor our conditions for the validity of Edgeworth expansions. Moreover, our proof for the first order also provides us with a result about Edgeworth expansions for smooth functions.

Schizophrenia is a debilitating neuropsychiatric disorder characterized by profound cortical network dysregulation, for which objective, clinically translatable EEG based biomarkers remain underdeveloped. Existing automated classification pipelines rely predominantly on static power spectral density features inherently blind to amplitude modulation dynamics and cross-frequency coupling, phenomena central to schizophrenia pathophysiology, while adopting epoch level cross validation strategies that introduce temporal data leakage, artificially inflate reported performance. This study introduces a mathematically principled diagnostic framework integrating the multi-order Wavelet Scattering Transform(WST), strict Leave One Subject Out (LOSO) cross-validation, and SHAP explainability for simultaneous EEG classification and biomarker discovery. Hierarchical WST coefficients capturing multi-scale amplitude modulation structure were extracted from resting state multichannel EEG. Subject-level ANOVA with Benjamini Hochberg false discovery rate correction identified significant biomarkers, with Random Forest and SVM classifiers evaluated under strict LOSO cross validation and subject-level majority voting. Second-order scattering coefficients encoding cross frequency coupling dominated the discriminative biomarker set, with gamma-band features most prevalent, demonstrating that temporal amplitude modulation constitutes the primary electrophysiological signature of schizophrenia. Electrode P3 was identified as the single most discriminative site. Under rigorous subject independent evaluation, the Random Forest achieved 90.48% accuracy (AUC = 0.9339; sensitivity = 95.56%). The proposed WST framework establishes a rigorous, interpretable standard for EEG-driven psychiatric biomarker discovery that can also be applicable in the detection of schizophrenia subtypes in the future.

We propose Emputation, a deep generative framework for learning imputation models. Emputation targets the extrapolation distribution of missing variables given observed variables, and training is guided by specific missingness assumptions that guarantee identification of the target distribution. The training objective, called the emputation risk, is an energy-score-based risk in which the identification assumption determines how observed entries are masked and which observations contribute to training. The resulting framework enables direct conditional sampling for multiple imputation. We show that the population minimizer of the emputation risk recovers the target extrapolation distribution under a broad class of identification assumptions, including several missing-not-at-random assumptions. Simulations show strong performance under both pointwise and distributional evaluation metrics, and an application to an Alzheimer's disease dataset demonstrates its practical value.

The $\beta$-model is popular for characterizing the commonly observed degree heterogeneity phenomenon in real-world networks. In this study, we develop a cycle counting approach to estimate $n$ node-specific parameters in the $\beta$-model for moderate or extremely sparse networks. Our proposed estimators, called Cycle Counting Ratio (CCR) Estimator, are based on the log-ratios of two network cycle counting statistics with explicit expressions and therefore easy to compute. We focus on conditions to guarantee statistical properties of the single estimator for each node. Under the very weak conditions that $\max_t \theta_t \to 0$ and $\theta_t \|\theta\|_1 \to \infty$, we show that the CCR estimator is consistent and achieves the minimax rate in terms of the mean squared error, which is the squared signal-to-noise ratio for $\hat{\beta}_t$ up to a constant factor. Here, $\hat{\beta}_t$ is the CCR estimator of the node-specific parameter $\beta_t$, $\theta_t = \exp(\beta_t)$ and $\theta=(\theta_1, \ldots, \theta_n)$. Even if the whole network density is close to the Erdős-Rényi lower bound $\log n/n$, the CCR estimator for the single parameter $\beta_t$ is still consistent as long as $\theta_t \|\theta\|_1 \to \infty$. To the best of our knowledge, this is the first time to derive the minimax rate and consistency result under such weak conditions. Under a slight stronger condition, we further establish its uniform consistency and asymptotic normality, whose asymptotic variance is $\theta_t \|\theta\|_1$. Numerical studies and an application to a sparse network data set demonstrate our theoretical findings.

We present msPCA: an open-source R package for sparse principal component analysis with multiple components. It implements an alternating maximization algorithm to generate a set of sparse loading vectors that collectively explain a large fraction of the variance in a dataset, while remaining non-redundant. The algorithm supports two definitions of non-redundancy: either orthogonality of the loading vectors or zero pairwise correlation between principal components (PCs). In the reported benchmarks, msPCA solves sparse PCA problems with thousands of features, achieving competitive runtimes while producing sparse components with controlled feasibility violations and a high fraction of variance explained.

We determine exactly what a kurtosis bound buys for one-sided tail control. For the class $\mathcal{C}(\kappa)$ of real random variables with mean $0$, variance $1$, and fourth moment at most $\kappa$, the skewness left free, we compute the worst-case tail probability $V_1(t,\kappa)=\sup_{X\in\mathcal{C}(\kappa)}\mathbb{P}(X\geq t)$ for every threshold $t>0$ and every $\kappa\geq 1$. The answer is a four-regime map: a Cantelli tongue $b(\kappa)\le t\le c(\kappa)$ on which the two-moment bound $1/(1+t^2)$ remains tight and the kurtosis constraint is worthless; a tail regime $t\geq c(\kappa)$ with the closed form $V_1=(\kappa-1)/((t^2-1)^2+\kappa-1)$; a plateau regime, present only for $\kappa\le 3/2$, on which the worst case freezes and the value does not depend on $t$; and a central regime described exactly by an explicit algebraic system, provably admitting no closed form in nested square roots. Beyond $c(\kappa)$ the one-sided and two-sided worst cases coincide: Cantelli's improvement over Chebyshev is annihilated by fourth-moment information. The minimal degree of a sum-of-squares proof of the tight bound is $2$ on the closed tongue and $4$ everywhere else, an exact phase diagram of proof degree. Every closed-form regime carries an explicit dual certificate and an explicit extremal distribution, re-verified on parameter grids by an independent checker in exact arithmetic. The closed forms invert to exact worst-case quantiles, sharpen a median-of-means constant, and give the exact per-direction tail available to degree-4 reasoning under certifiable kurtosis. We found the map through an AI-guided search around the certifying pipeline, LemmaForge, which is validated on classical benchmarks, independently reproduces the symmetric-slice bound of Zelen (1954), and recovers the $2\sqrt{3}-3$ constant of He, Zhang, and Zhang (2010) at $t=0$.

We study minimax goodness-of-fit testing for uniformity from $n$ multinomial observations over $N$ categories against $\ell_p$ departures of size $\epsilon_n$. Writing $u_n:=\epsilon_n^2 n\,N^{3/2-2/p}/\sqrt{2}$ for the associated signal-to-noise ratio, we focus on the intermediate regime $N=o(n^2)$ with $u_n\to u^*\in(0,\infty)$, in which the minimax risk converges to a nontrivial constant. In the Poissonized version of the problem this constant equals $2\Phi(-u^*/2)$ \cite{Kipnis2025minimax}, yielding an upper bound on the multinomial minimax risk. Here we prove the matching lower bound. The key step is a conditional central limit theorem for weighted sums under a Poisson mixture prior, conditioned on the total count. Together with the upper bound in \cite{Kipnis2025minimax}, this gives an exact sharp-constant characterization of the multinomial minimax risk in the intermediate regime.

Diabetic retinopathy (DR) is a local retinal lesion process and a visible manifestation of systemic microvascular injury. Modern retinal AI can grade images accurately, but often leaves unanswered how local lesion evidence, retinal vascular structure, and systemic disease pathways are connected. This paper introduces Causal-RetiGraph, a compact biomedical informatics framework that links retinal graph phenotypes with NHANES-anchored pathway modelling. The retinal-image fold constructs an interpretable $X1234$ phenotype from vessel maps, lesion evidence, image embeddings, and AutoMorph biomarkers through spatial $X_{12}$ and Jacobian $X_{34}$ branches. The NHANES fold models systemic exposures, covariates, a same-subject retinal mediator family $R^*$, and downstream outcome families. $X1234$ is used for retinal support and pathway prioritisation, while $R^*$ is used for participant-level pathway summaries. On the retinal fold, $X1234$ achieves 0.9055 binary DR accuracy and 0.9711 AUROC, with graded DR QWK of 0.8312. The results show that lesion and biomarker streams improve contextual retinal representation under scarce and imbalanced data. In NHANES, HbA1c, urine albumin, pulse pressure, fasting glucose, and systolic blood pressure are the strongest binary DR anchors. Participant-level pathway analysis identifies glycaemic–renal and glycaemic–haemodynamic pathways as the clearest mediator-style signals. These results suggest that retinal graph phenotypes can help prioritise systemic pathways in DR while preserving the distinction between image-derived support and same-subject mediation.

Scientists often seek to draw causal inferences from structured data that is not independently and identically distributed, such as spatial data, network data, or molecular data. We develop geometric causal models (GCMs), a framework for causal inference from dependent data that exploits underlying symmetries of the data generating process. For example, in spatial data, we consider processes that are symmetric under translations, or in graph data, symmetric under permutations of the nodes. We show how symmetries, formalized with group theory, can enable causal identification and estimation. We deploy ergodic theory for amenable groups to establish identification, and combine geometric deep learning with scalable Bayesian inference for estimation. We recover i.i.d. causal models and do-calculus when the data is a sequence and the symmetry is permutation equivariance, and find novel types of causal models when we use alternate structures and symmetries. As an example, we construct a causal model that satisfies the symmetries of DNA. This GCM enables new estimators for the effects of genetic variation, combining deep functional genomics models to describe outcomes and DNA language models to describe propensities. We illustrate on semisynthetic data.

We focus on regression settings where the fairness goal is to equalize average predictions across values of a sensitive attribute, a criterion known as mean demographic parity. Directly optimizing this criterion is difficult because it depends on a conditional mean that is unknown and changes during training. Common dependence penalties and adversarial methods do not estimate this conditional mean; instead, they push predictions toward full independence. This stronger constraint can reduce accuracy even when average predictions are already equal. Existing conditional-mean methods are limited to linear predictors or low-dimensional sensitive attributes. We enable direct optimization of mean demographic parity using DPVar, a fairness measure defined as the variance of the conditional mean prediction. Because the conditional mean must be estimated as the predictor changes, optimizing DPVar leads to a functional bilevel problem. We develop two solvers: FBO, which uses a closed-form hypergradient, and an iterative-differentiation (ITD) solver that differentiates through updates of the conditional-mean estimator. Unlike previous conditional-mean methods, our approach applies to nonlinear predictors and high-dimensional continuous sensitive attributes. Across a semi-synthetic benchmark built from 21 tabular regression datasets and Communities & Crime data, FBO and ITD recover competitive or better accuracy-DPVar trade-offs than existing methods.

Linear regression is one of the core statistical tools used for analysis of data. In the era of statistical learning, linear regression has been expanded into two directions. The first is regularisation, where penalties are added to the loss function to obtain more stable or sparse solutions. The second direction is basis expansion, such as with spline or kernel functions, where the linearity assumption is dropped. In practice, empirical researchers often collect multiple outcome variables. Regression models, either linear, regularized, or expanded, can be fitted to each of these outcome variables, but such an approach does not take into account the associations among the response variables. Reduced rank regression is a multivariate regression tool that takes into account and models the association among the response variables. In this paper, we develop and test a B-spline basis expansion for reduced rank regression, a regression model for multiple outcome variables where we penalize the coefficients to obtain a smooth fit, as in P-splines. This approach is useful for the analysis of data sets with multiple outcomes and a relatively small number of predictors. A block-relaxation algorithm is developed for parameter estimation and we discuss ways for penalty parameter tuning. We develop visualizations tools for model interpretation based on biplot methodology for “all predictor - all response” relationships and partial dependence plots for interpreting “single predictor - single response” relationships. With show several experimental gauges and we analyze three empirical data sets. We conclude this paper with a discussion.

Monetary risk measures have gained popularity for expressing decision-makers' risk aversion. Value-at-Risk (VaR) and Conditional-Value-at-Risk (CVaR), in particular, are used commonly for this purpose. This paper proposes new efficient algorithms to compute these risk measures for a discrete random variable in expected linear time with respect to the size of its domain. First, we propose a QuickVaR algorithm that computes the VaR of a discrete random variable. Then, we leverage QuickVaR to propose QuickDivergence, an algorithm for computing a class of $\varphi$-divergence risk measures, including the popular CVaR risk measure. The QuickVaR algorithm adapts the well-known Quickselect algorithm, while QuickDivergence builds on polymatroid optimization algorithms. Numerical results show that our new algorithms offer an order-of-magnitude speedup for large domains, and a library implementation of the algorithms is available at https://github.com/RiskAverseRL/RiskMeasures.jl.

This paper investigates the problem of statistical inference for a mixture distribution consisting of a discrete and a continuous component, with a particular focus on the class of rational-infinitely divisible distributions. We consider non-parametric estimation of both components of the mixture as well as the quasi-Lévy measure, assuming that the mixture belongs to the class of rational-infinitely divisible distributions. We propose an estimation framework based on band-limited kernels, which are the functions characterized by compactly supported Fourier transform. Under mild assumptions, the proposed estimators are theoretically shown to achieve polynomial (and in some cases even almost parametric) convergence rates. Finally, we demonstrate the numerical performance of the algorithm on simulated examples.

We study post-learning inference for structural properties of data-dependent combinatorial optimizers. The target is whether an oracle optimizer, rather than a latent parameter or smooth functional, belongs to a prescribed class, such as a category-mix, inventory, or resource-feasibility class. We focus on a high-dimensional contextual multinomial logit model with sequentially adaptive data collection, where the parameter-to-optimizer map is discontinuous and the policy induces temporal dependence. We propose a novel perturbation test based on a nonsmooth max-difference revenue statistic comparing the best null assortment with the best alternative assortment. The test perturbs the estimated terminal revenue surface on the selected support: random unit directions capture directional uncertainty, while the minimal perturbation radius captures magnitude uncertainty and yields a p-value. This localizes inference near the null–alternative boundary and avoids uniform error control over the full candidate class. The data are collected by an \(\ell_1\)-penalized online likelihood policy that performs variable selection while controlling regret. Using a new anti-concentration argument for Gaussian maxima differences and martingale Gaussian coupling, we establish uniform estimation rates, effective support recovery, and asymptotic validity of the proposed p-value under adaptive assortment selection. We prove asymptotic size control and power consistency under a localized signal condition.

Accurately estimating the unknown target label distribution is the critical first step for adapting to label shift. This task, widely known as quantification or class prevalence estimation, has recently seen significant advances through continuous KDE-based methods which model the density of multiclass classifier posteriors. Posterior vectors might be regarded as compositional data, since they lie on the probability simplex. However, existing KDE-based quantifiers typically rely on Euclidean Gaussian kernels, which ignore simplex geometry and incorrectly assign probability mass outside its boundaries. We introduce a geometry-aware KDE model for multiclass quantification based on log-ratio representations and Aitchison geometry, together with a shrinkage regularization that improves robustness near the simplex boundary. Combined with a maximum-likelihood interpretation of KDE-based quantification, we derive both point-estimation and Bayesian inference procedures for class prevalences. Experiments on 42 datasets across tabular, text, and image domains show that the proposed method is competitive with state-of-the-art quantifiers, often improving over standard KDE-based baselines, while also yielding strong results among Bayesian quantification methods.

While Bayesian methods are increasingly used in clinical research, confusion persists as to whether Bayesian designs are affected by repeated interim analyses, and how such effects should be evaluated. We aimed to clarify this question by evaluating both frequentist and Bayesian operating characteristics of group-sequential trials. We conducted simulation studies with normally distributed outcomes, examining designs with repeated analyses, with and without futility stopping rules and multiplicity adjustments. We show that repeated interim analyses alter both frequentist and Bayesian operating characteristics in group-sequential trials, regardless of the inferential framework adopted. Without proper adjustment for multiplicity, the Type I error rate increases with the number of analyses. Bayesian operating characteristics such as the risk of erroneous conclusions and the informative value of an efficacy conclusion are meaningful alternatives to classical frequentist metrics, but they are sensitive to prior divergence between stakeholders, and this sensitivity grows with the number of analyses. Even under full prior agreement, the informative value of an efficacy conclusion is reduced. While it is possible to control frequentist operating characteristics of Bayesian trials with appropriate methods, it is not possible to guarantee such control for Bayesian operating characteristics because they are prior-dependent and different stakeholders may adopt different priors. Contrary to claims that Bayesian inference is immune to multiplicity, our results show that Bayesian clinical trials are no less affected by repeated analyses than frequentist ones, regardless of the framework used to evaluate them.

Evaluating the causal effect of a treatment on an outcome is a central objective in causal inference. While the average causal effect summarizes the mean impact of treatment, the central moments of the individual causal effect (ICE) characterize the shape of the ICE distribution, thereby revealing the extent and structure of treatment effect heterogeneity across individuals. This paper investigates the identification and bounding of the central moments of the ICE using only the marginal central moments of each potential outcome (PO). Compared with existing approaches that require knowledge of the full marginal distributions of the POs, marginal moment information is often substantially easier to obtain in empirical applications. Finally, we illustrate the practical relevance of our results through two empirical case studies.

Given a nonnegative matrix $X$, a factorization rank $r$ and {a positive integer $p$}, entrywise power matrix factorization (EPMF) looks for a low-rank matrix $X_r$ such that $X = |X_r|^{\circ p}$ (exact case) or $X \approx |X_r|^{\circ p}$ (approximate case), where $(\cdot)^{\circ p}$ denotes the componentwise exponent. EPMF includes the modulus model ($p=1$) and componentwise square factorization ($p=2$) as special cases, the latter being closely related to the square root rank. We analyze the computational complexity of the exact decision problem and the Frobenius-norm approximation problem, and establish a complete complexity landscape. In the exact case, we show that EPMF is equivalent to the combinatorial problem of flipping the signs of the entries of a given matrix $X$ to obtain a rank-$r$ matrix, which we refer to as the low-rank matrix signing (LRMS) problem. We first show that LRMS, and hence exact EPMF, is strongly NP-hard, improving a weak NP-hardness result for the square-root-rank (Math. Prog., 2015). We then show that LRMS can be solved in polynomial time when $r$ is fixed. Moreover, when the rank $r$ is part of the input, we show that for generic matrices the algorithm is fixed-parameter tractable (FPT) in the parameter $r$; in fact, the running time is fixed-parameter linear in the number of entries of the input matrix. In the approximate case using the Frobenius norm as an error measure, we show that EPMF is NP-hard, already when $r=2$, the smallest nontrivial case.

Motivated by real-world scenarios where malicious entities tamper with existing networks, we define a model where an adversary seeks to hide a set of corrupted vertices inside a graph $G^*$. To this end, the adversary can add edges between the corrupted vertices, as well as edges between the corrupted vertices and $G^*$, and its power is then measured by the size of the neighborhood of the corrupted vertices in $G^*$. Our goal is to design an active learning algorithm that efficiently finds the subset of corrupted vertices using a small number of label queries. We devise an efficient algorithm that approximately recovers the corrupted vertices with a query complexity that depends polynomially on both the power of the adversary and the vertex expansion of $G^*$, a fundamental measure of graph connectivity. At the heart of this result is a polynomial-time algorithm, obtained by carefully adapting sum-of-squares algorithms for approximating minimum expansion, that finds a set with small vertex expansion subject to cardinality constraints. To the best of our knowledge, this is the first time that the vertex expansion is shown to play a key role in determining the query complexity of active learning algorithms robust to structural adversarial attacks.

We present a Bayesian P-spline method for estimating the frequency-dependent cross-spectral density matrix of stationary multivariate time series. The inverse spectral matrix is parametrised through its frequency-varying Cholesky decomposition, which guarantees Hermitian positive definiteness at every frequency. Each real log-diagonal entry and each real and imaginary off-diagonal entry is given an independent penalised B-spline prior that controls smoothness. Inference uses a blocked, coarse-grained Whittle likelihood with safe-Bayes $\eta$-tempering to stabilise posterior calibration, sampled by the No-U-Turn Sampler from a variational initialisation. On synthetic VAR(2) benchmarks with known ground truth, the method recovers both diagonal and cross-spectral structure, attains near-nominal credible-interval coverage, and achieves a relative integrated squared (Frobenius) error (RISE) that decreases with sample size. We then apply the method to publicly released simulated LISA time-delay interferometry (TDI) data in two noise configurations. In the idealised symmetric case, the full multivariate model and a reduced model that assumes a diagonal AET noise covariance agree to within $\sim10^{-3}$ in RISE. Under realistic noise that is asymmetric across the six Movable Optical Sub-Assemblies (MOSAs), the AET-diagonal assumption fails by more than an order of magnitude in RISE ($\sim\!3.3\!\times\!10^{-2}$ versus $\sim\!10^{-3}$), whereas the full multivariate model recovers the cross-spectral structure.

We study a sequential learning problem for stable matchings in two-sided markets where preferences on both sides are initially unknown. We focus on a centralized setting where an algorithm matches agents at each time step and receives noisy rewards that reflect the preferences of the matched agents, following a semi-bandit feedback structure. We adopt a pure exploration perspective, aiming to efficiently identify the optimal stable matching with high probability. Our work extends prior results by handling two-sided uncertainty and by exploiting partial preference information. A central ingredient is the notion of pervasive stable matching, which enables the identification of optimal stable matchings under partial preferences. We propose elimination-based algorithms whose stopping criteria exploit the structure of the learned partial preferences, and provide a refined sample-complexity analysis. Beyond pure exploration, we extend our approach to regret minimization and establish regret bounds with respect to the optimal stable matching that avoid dependence on the minimum reward gap $\Delta_{\min}$.

The Functional Autoregressive Model (FAR) generalizes the multivariate AR(1) model in Time Series Analysis to functional data. It serves as a historical foundational point in the study of functional time series and remains a fundamental and widely used model for dependent functional data. The process-the observed data-generated by the FAR model forms a Hilbert-valued Markov chain. This paper investigates the non-asymptotic prediction mean square error and derives a lower bound. This lower bound is established in the specific context of non-i.i.d. data and depends on the mixed smoothness of the functional time series and of the unknown correlation operator driving the FAR model. Instead of the standard functional PCA regularization, a ridge-type estimator is proposed, which avoids the preliminary estimation of the spectrum of the covariance sequence associated with the process. A non-asymptotic upper bound is derived for this estimate, which matches the lower bound up to multiplicative constants. Furthermore, a detailed study of the estimate's bias reveals connections between functional smoothness parameters and regularly/rapidly varying functions, which are common in extreme value theory. Simulation results corroborate the theoretical main theorems.

Tabular Foundation Models (TFMs) have demonstrated strong empirical performance as black-box inference engines through in-context learning. However, their use in transfer learning is limited by two obstacles: strict context-size constraints and sensitivity to distribution shifts between source and target tasks. Directly pooling heterogeneous source data can therefore lead to negative transfer. To address these challenges, we propose Context-Constrained Transfer Learning via ANchoring and DIstillation (TL-ANDI), a posterior-aware distillation framework for TFMs. TL-ANDI constructs a compact source context by solving a budget-constrained optimal transport problem whose cost jointly measures target covariate coverage and posterior compatibility. The selected anchor samples are then equipped with locally distilled labels and combined with a residual calibration step using target data.

In hybrid uncertainty quantification, evaluating how aleatory sensitivities vary under epistemic uncertainty, referred to as conditional Sobol' indices, is typically hindered by the computationally expensive double-loop procedure. Classical Polynomial Chaos Expansion (PCE) provides efficient access to global sensitivity measures but cannot directly resolve sensitivity variation across the epistemic space without repeated surrogate reconstruction. This study proposes a unified Bayesian framework that extracts continuous conditional Sobol' fields from a single global PCE representation. By exploiting the tensor-product structure of orthogonal polynomial bases in an augmented stochastic space, the global expansion is analytically decomposed into epistemic-dependent coefficient fields, enabling a closed-form variance decomposition. As a result, conditional Sobol' indices can be computed through a purely algebraic post-processing step without additional model evaluations or retraining. In addition, a Reversible Jump Markov Chain Monte Carlo (RJMCMC) scheme is incorporated to perform adaptive basis selection and trans-dimensional inference, while simultaneously providing Bayesian credible intervals for the conditional sensitivity measures. Numerical experiments on a high-dimensional groundwater flow model demonstrate that the proposed method significantly reduces computational cost while maintaining smooth sensitivity fields and statistically consistent uncertainty quantification across the epistemic domain.

Post-hoc conditioning of pretrained diffusion models can be addressed using Sequential Monte Carlo (SMC) methods. By evolving an interacting particle system, SMC-guided diffusion samplers combine unconditional reverse-diffusion dynamics with sequential reweighting to approximate conditional distributions. Nevertheless, even in the infinite-particle limit, the implemented sampler may differ from the ideal conditional target because of errors in the diffusion model, its numerical implementation, and the guidance mechanism. We characterize how these local errors propagate through forward-smoothing kernels, which jointly account for the reverse dynamics and the remaining conditioning information. This yields non-asymptotic error bounds that capture both finite-particle fluctuations and approximation errors arising from initialization, numerical integration, score approximation, and potential design. In doing so, we extend stability guarantees for diffusion models to the conditional setting. Finally, we apply our framework to several state-of-the-art SMC-guided diffusion algorithms, providing a unified theoretical perspective on their approximation mechanisms and sources of error.

Score-based Generative Models (SGMs) have achieved impressive performance in data generation across a wide range of applications. While the statistical properties of their sampling procedures are increasingly well understood, the optimization dynamics underlying their training remain less explored. SGMs are typically trained by minimizing a weighted denoising score-matching objective, yet optimization guarantees with stochastic gradients remain limited. In this work, we study Stochastic Gradient Descent (SGD) for SGMs, contributing results in two complementary regimes. For general score parameterizations, we derive a non-convex analysis of SGD for the weighted denoising score-matching objective, making explicit how the resulting optimization bound depends on the loss weighting and time-sampling distribution. We then consider overparameterized two-layer ReLU networks and develop a Neural Tangent Kernel analysis tailored to diffusion training with stochastic gradients, yielding score-approximation error bounds along the SGD trajectory. Our analysis quantifies the role of the reweighting factor in these bounds, providing a theoretical characterization of weighting choices used in practice.

We study on-policy distillation (OPD) for agentic tasks, where an LLM agent interacts with an environment over multiple turns and a student imitates a teacher over these multi-turn interaction histories. Fully online OPD is costly because each update requires fresh student rollouts through the environment and teacher queries at visited histories. We propose Replayed-Prefix On-Policy Distillation (ReOPD), an off-environment alternative that reuses pre-collected teacher trajectories as replayed prefixes: the student acts at selected steps, while the teacher provides dense per-step supervision without executing new environment interactions. We show that multi-turn OPD introduces a prefix trap: making histories more student-on-policy improves relevance to the student, but can query the teacher on histories where its target is unreliable. This creates a two-sided distribution shift between student occupancy and teacher reliability. ReOPD addresses this by treating multi-turn OPD as a reliability-aware prefix distribution design and implements it with a simple step-decaying sampling schedule that emphasizes early, lower-shift prefixes. Across mathematical reasoning with Python and search environments over multiple teacher and student model scales, ReOPD preserves or improves OPD-level accuracy, uses zero tool calls during student training, and is at least 4$\times$ faster per rollout than OPD. ReOPD therefore turns expensive agent-environment interaction into a reusable offline resource, enabling scalable distillation across tools, tasks, and environments.

Higher-order influence functions, introduced in a series of articles (Robins et al., 2008, 2009a; van der Vaart, 2014; Robins et al., 2016, 2023; Liu et al., 2017), are a unified framework for constructing rate-optimal point estimates of a class of statistical functionals under various complexity-reducing assumptions on the posited statistical model that generates the observed data. Although higher-order (influence functions) estimators are theoretically appealing, they have very limited practical uptake compared to their first-order counterparts. The original higher-order estimators proposed in Robins et al. (2008) and Robins et al. (2017) involve nonparametric density estimation of multi-dimensional covariates, a highly nontrivial statistical and computational problem on its own. The density estimator is, in turn, used in the evaluation of the inverse population Gram matrix $\Omega$ of a set of $k$-dimensional basis transformations of covariates. There, $k$ is allowed to be as large as $o (n^2)$. To partially address this potential shortcoming, Liu et al. (2017) restrict $k$ to $o (n)$ and instead estimate $\Omega$ directly using the inverse sample Gram matrix estimator, but computed from an independent sample often obtained by sample-splitting. Liu et al. (2017) refer to this alternative estimator as the empirical higher-order estimator. Although the empirical higher-order estimator bypasses density estimation, it suffers from numerical instability due to inverting a large-dimensional sample Gram matrix. In this article, for a class of bilinear forms/functionals that often appear in substantive fields, we propose a new stabilized higher-order estimator without sample splitting, which exhibits more stable finite-sample performance compared to the empirical higher-order estimator. We also prove that this new class of higher-order estimators enjoys similar statistical guarantees.

Reconstructing population dynamics is a central problem in the physical and data sciences. Often, the dynamics are modeled as a Wasserstein gradient flow (WGF): a curve of distributions driven by an energy functional. Though there are multiple mathematical characterizations of a WGF, the dominant algorithmic approach relies on the Jordan–Kinderlehrer–Otto (JKO) scheme. JKO-based methods are inflexible to time discretisation and require solving costly optimal transport problems. We take a residual approach, enforcing the continuity equations via a non-negative loss function whose minimum is the WGF. Combined with a data-fitting divergence, this gives a single global objective. This perspective unifies several existing methods and leads to a new particle-based method, stitching, that is simulation-free and robust to large gaps between observations. We demonstrate that the stitching method achieves state-of-the-art performance across trajectory inference benchmarks. For code see github.com/BasisResearch/wasserstein-residuals.