Papers

Machine learning methods estimate flexible relationships between confounders, treatments and outcomes, making unconfoundedness assumptions in causal analysis more palatable. How successful are these approaches in recovering ground truth baselines? We analyze new data covering an experimental rollout of a software product feature for one group of users and a voluntary opt-in to the same feature for another similar group. Recovering causal effects from nonrandomized data is feasible–but only with careful modeling choices. We offer best practices for more credible treatment effect estimation in modern, high-dimensional datasets and evaluate which tests provide meaningful signals about the reliability of estimated effects.

Global sensitivity analysis of complex numerical simulators is often limited by the small number of model evaluations that can be afforded. In such settings, surrogate models built from a limited set of simulations can substantially reduce the computational burden, provided that the design of computer experiments is enriched efficiently. In this context, we propose an active learning approach that, for a fixed evaluation budget, targets the most informative regions of the input space to improve sensitivity analysis accuracy. More specifically, our method builds on recent advances in active learning for sensitivity analysis (Sobol' indices and derivative-based global sensitivity measures, DGSM) that exploit derivatives obtained from a Gaussian process (GP) surrogate. By leveraging the joint posterior distribution of the GP gradient, we develop acquisition functions that better account for correlations between partial derivatives and their impact on the response surface, leading to a more comprehensive and robust methodology than existing DGSM-oriented criteria. The proposed approach is first compared to state-of-the-art methods on standard benchmark functions, and is then applied to a real environmental model of pesticide transfers.

We study nonasymptotic (finite-sample) confidence intervals for treatment effects in randomized experiments. In the existing literature, the effective sample sizes of nonasymptotic confidence intervals tend to be looser than the corresponding central-limit-theorem-based confidence intervals by a factor depending on the square root of the propensity score. We show that this performance gap can be closed, designing nonasymptotic confidence intervals that have the same effective sample size as their asymptotic counterparts. Our approach involves systematic exploitation of negative dependence or variance adaptivity (or both). We also show that the nonasymptotic rates that we achieve are unimprovable in an information-theoretic sense.

Explicit modelling of between-study heterogeneity is essential in network meta-analysis (NMA) to ensure valid inference and avoid overstating precision. While the additive random-effects (RE) model is the conventional approach, the multiplicative-effect (ME) model remains underexplored. The ME model inflates within-study variances by a common factor estimated via weighted least squares, yielding identical point estimates to a fixed-effect model while inflating confidence intervals. We empirically compared RE and ME models across NMAs of two-arm studies with significant heterogeneity from the nmadb database, assessing model fit using the Akaike Information Criterion. The ME model often provided comparable or better fit to the RE model. Case studies further revealed that RE models are sensitive to extreme and imprecise observations, whereas ME models assign less weight to such observations and hence exhibit greater robustness to publication bias. Our results suggest that the ME model warrant consideration alongside conventional RE model in NMA practice.

We consider the problem of learning the network of mutual excitations (i.e., the dependency graph) in a non-stationary, multivariate Hawkes process. We consider a general setting where baseline rates at each node are time-varying and delay kernels are not shift-invariant. Our main results show that if the dependency graph of an $n$-variate Hawkes process is sparse (i.e., it has a maximum degree that is bounded with respect to $n$), our algorithm accurately reconstructs it from data after observing the Hawkes process for $T = \mathrm{polylog}(n)$ time, with high probability. Our algorithm is computationally efficient, and provably succeeds in learning dependencies even if only a subset of time series are observed and event times are not precisely known.

Algorithmic stability is a central concept in statistics and learning theory that measures how sensitive an algorithm's output is to small changes in the training data. Stability plays a crucial role in understanding generalization, robustness, and replicability, and a variety of stability notions have been proposed in different learning settings. However, while stability entails desirable properties, it is typically not sufficient on its own for statistical learning – and indeed, it may be at odds with accuracy, since an algorithm that always outputs a constant function is perfectly stable but statistically meaningless. Thus, it is essential to understand the potential statistical cost of stability. In this work, we address this question by adopting a statistical decision-theoretic perspective, treating stability as a constraint in estimation. Focusing on two representative notions-worst-case stability and average-case stability-we first establish general lower bounds on the achievable estimation accuracy under each type of stability constraint. We then develop optimal stable estimators for four canonical estimation problems, including several mean estimation and regression settings. Together, these results characterize the optimal trade-offs between stability and accuracy across these tasks. Our findings formalize the intuition that average-case stability imposes a qualitatively weaker restriction than worst-case stability, and they further reveal that the gap between these two can vary substantially across different estimation problems.

Dense Associative Memory (DAM) models generalize the classical Hopfield model by incorporating n-body or exponential interactions that greatly enhance storage capacity. While the criticality of DAM models has been largely investigated, mainly within a statistical equilibrium picture, little attention has been devoted to the temporal self-organizing behavior induced by learning. In this work, we investigate the behavior of a stochastic exponential DAM (SEDAM) model through the lens of Temporal Complexity (TC), a framework that characterizes complex systems by intermittent transition events between order and disorder and by scale-free temporal statistics. Transition events associated with birth-death of neural avalanche structures are exploited for the TC analyses and compared with analogous transition events based on coincidence structures. We systematically explore how TC indicators depend on control parameters, i.e., noise intensity and memory load. Our results reveal that the SEDAM model exhibits regimes of complex intermittency characterized by nontrivial temporal correlations and scale-free behavior, indicating the spontaneous emergence of self-organizing dynamics. Notably, such regimes arise over finite intervals of noise intensity rather than at a single critical point, consistent with the concept of extended criticality. Further, the noise intensity range needed to reach the critical region, where self-organizing behavior emerges, slightly decreases as the memory load increases. This study highlights the relevance of TC as a complementary framework for understanding learning and information processing in artificial and biological neural systems, revealing the link between the memory load and the self-organizing capacity of the network.

This article is an exposition on some recent theoretical advances in learning latent structured models, with a primary focus on the fundamental roles that optimal transport distances play in the statistical theory. We aim at what may be the most critical and novel ingredient in this theory: the motivation, formulation, derivation and ramification of inverse bounds, a rich collection of structural inequalities for latent structured models which connect the space of distributions of unobserved structures of interest to the space of distributions for observed data. This theory is illustrated on classical mixture models, as well as the more modern hierarchical models that have been developed in Bayesian statistics, machine learning and related fields.

We introduce a smooth variant of the SCAD thresholding rule for wavelet denoising by replacing its piecewise linear transition with a raised cosine. The resulting shrinkage function is odd, continuous on R, and continuously differentiable away from the main threshold, yet retains the hallmark SCAD properties of sparsity for small coefficients and near unbiasedness for large ones. This smoothness places the rule within the continuous thresholding class for which Stein's unbiased risk estimate is valid. As a result, unbiased risk computation, stable data-driven threshold selection, and the asymptotic theory of Kudryavtsev and Shestakov apply. A corresponding nonconvex prior is obtained whose posterior mode coincides with the estimator, yielding a transparent Bayesian interpretation. We give an explicit SURE risk expression, discuss the oracle scale of the optimal threshold, and describe both global and level-dependent adaptive versions. The smooth SCAD rule therefore offers a tractable refinement of SCAD, combining low bias, exact sparsity, and analytical convenience in a single wavelet shrinkage procedure.

Diffusion models now generate high-quality, diverse samples, with an increasing focus on more powerful models. Although ensembling is a well-known way to improve supervised models, its application to unconditional score-based diffusion models remains largely unexplored. In this work we investigate whether it provides tangible benefits for generative modelling. We find that while ensembling the scores generally improves the score-matching loss and model likelihood, it fails to consistently enhance perceptual quality metrics such as FID on image datasets. We confirm this observation across a breadth of aggregation rules using Deep Ensembles, Monte Carlo Dropout, on CIFAR-10 and FFHQ. We attempt to explain this discrepancy by investigating possible explanations, such as the link between score estimation and image quality. We also look into tabular data through random forests, and find that one aggregation strategy outperforms the others. Finally, we provide theoretical insights into the summing of score models, which shed light not only on ensembling but also on several model composition techniques (e.g. guidance).

Gaussian Random Fields (GRFs) with Matérn covariance functions have emerged as a powerful framework for modeling spatial processes due to their flexibility in capturing different features of the spatial field. However, the smoothness parameter is challenging to estimate using maximum likelihood estimation (MLE), which involves evaluating the likelihood based on the full covariance matrix of the GRF, due to numerical instability. Moreover, MLE remains computationally prohibitive for large spatial datasets. To address this challenge, we propose the Fisher-BackTracking (Fisher-BT) method, which integrates the Fisher scoring algorithm with a backtracking line search strategy and adopts a series approximation for the modified Bessel function. This method enables an efficient MLE estimation for spatial datasets using the ExaGeoStat high-performance computing framework. Our proposed method not only reduces the number of iterations and accelerates convergence compared to derivative-free optimization methods but also improves the numerical stability of the smoothness parameter estimation. Through simulations and real-data analysis using a soil moisture dataset covering the Mississippi River Basin, we show that the proposed Fisher-BT method achieves accuracy comparable to existing approaches while significantly outperforming derivative-free algorithms such as BOBYQA and Nelder-Mead in terms of computational efficiency and numerical stability.

This work presents the first systematic development of Stein's method for matrix distributions. We establish the basic essential ingredients of Stein's method for matrix normal approximation: we derive an extended-generator-based Stein identity from a matrix Ornstein-Uhlenbeck diffusion with two-sided scales, provide an explicit semigroup representation for the solution of the Stein equation, and obtain regularity estimates for the solution. The new methodology is demonstrated in three examples: (i) smooth Wasserstein distance bounds to quantify the matrix central limit theorem (a didactic example), (ii) a Wasserstein distance bound for the matrix normal approximation of the centered matrix $T$ distribution, and (iii) a Stein's method-of-moments approach to estimating the row and column covariance factors of the matrix normal, yielding a flexible class of weighted flip-flop Stein estimators that generalize Dutilleul's classical flip-flop algorithm and naturally accommodate row/column importance weights, systematic missingness, and projection onto structured covariance families. The latter two examples are intrinsically matrix-valued and cannot be treated using naive vectorization.

Robustness under perturbation and contamination is a prominent issue in statistical learning. We address the robust nonlinear regression based on the so-called interval conditional value-at-risk (In-CVaR), which is introduced to enhance robustness by trimming extreme losses. While recent literature shows that the In-CVaR based statistical learning exhibits superior robustness performance than classical robust regression models, its theoretical robustness analysis for nonlinear regression remains largely unexplored. We rigorously quantify robustness under contamination, with a unified study of distributional breakdown point for a broad class of regression models, including linear, piecewise affine and neural network models with $\ell_1$, $\ell_2$ and Huber losses. Moreover, we analyze the qualitative robustness of the In-CVaR based estimator under perturbation. We show that under several minor assumptions, the In-CVaR based estimator is qualitatively robust in terms of the Prokhorov metric if and only if the largest portion of losses is trimmed. Overall, this study analyzes robustness properties of In-CVaR based nonlinear regression models under both perturbation and contamination, which illustrates the advantages of In-CVaR risk measure over conditional value-at-risk and expectation for robust regression in both theory and numerical experiments.

We derive explicitly the e-values with optimal (relative) growth rate in the worst case for testing the mean of a bounded random variable, thereby providing the first application of the (RE)GROW quality criteria beyond the assumption of mutually absolutely continuous hypotheses for e-values originally proposed by Grünwald et al. (2024). For both criteria, we explicitly characterise the alternatives that are most difficult to test against and show that they admit a meaningful interpretation. We give two important examples in which REGROW provides a powerful quality criterion to choose optimal e-variables whereas GROW leads to trivial solutions.

Hidden variable graphical models can sometimes imply constraints on the observable distribution that are more complex than simple conditional independence relations. These observable constraints can falsify assumptions of the model that would otherwise be untestable due to the unobserved variables and can be used to constrain estimation procedures to improve statistical efficiency. Knowing the complete set of observable constraints is thus ideal, but this can be difficult to determine in many settings. In models with categorical observed variables and a joint distribution that is completely characterized by linear relations to the unobservable response function variables, we develop a systematic method for deriving the complete set of observable constraints. We illustrate the method in several new settings, including ones that imply both inequality and equality constraints.

Time series often exhibit non-ergodic behaviour that complicates forecasting and inference. This article proposes a likelihood-based approach for estimating ergodicity transformations that addresses such challenges. The method is broadly compatible with standard models, including Gaussian processes, ARMA, and GARCH. A detailed simulation study using geometric and arithmetic Brownian motion demonstrates the ability of the approach to recover known ergodicity transformations. A further case study on the large macroeconomic database FRED-QD shows that incorporating ergodicity transformations can provide meaningful improvements over conventional transformations or naive specifications in applied work.

We define two minimum distance estimators for dependent data by minimizing some approximated Maximum Mean Discrepancy distances between the true empirical distribution of observations and their assumed (parametric) model distribution. When the latter one is intractable, it is approximated by simulation, allowing to accommodate most dynamic processes with latent variables. We derive the non-asymptotic and the large sample properties of our estimators in the context of absolutely regular/beta-mixing random elements. Our simulation experiments illustrate the robustness of our procedures to model misspecification, particularly in comparison with alternative standard estimation methods.

Qualitative Comparative Analysis (QCA) requires researchers to choose calibration and dichotomization thresholds, and these choices can substantially affect truth tables, minimization, and resulting solution formulas. Despite this dependency, threshold sensitivity is often examined only in an ad hoc manner because repeated analyses are time-intensive and error-prone. We present ThSQCA, an R package that automates threshold-sweep analyses by treating thresholds as explicit analytical variables. It provides four sweep functions (otSweep, ctSweepS, ctSweepM, dtSweep) to explore outcome thresholds, single-condition thresholds, multi-condition threshold grids, and joint outcome-condition threshold spaces, respectively. ThSQCA integrates with the established CRAN package QCA for truth table construction and Boolean minimization, while returning structured S3 objects with consistent print/summary methods and optional detailed results. The package also supports automated Markdown report generation and configuration-chart output to facilitate reproducible documentation of cross-threshold results.

Here, we introduce a stochastic partial differential equation (SPDE) formulation driven by temporally correlated noise to describe light propagation beyond the standard Markov approximation. By representing the squared refractive index fluctuations as a random field with explicit long-range temporal correlations, we demonstrate that the propagation dynamics map exactly onto the hyperbolic Anderson model. This rigorous mapping enables the derivation of new quantitative scaling relations that connect the environment's non-Markovian memory effects to the statistical properties of the emergent light field. We experimentally validate these analytical predictions in an outdoor atmospheric environment, confirming the memory-dependent statistical signatures of the propagated light. Our results establish a precise physical foundation for understanding memory-driven wave phenomena, providing crucial insights for free-space optical communication, remote sensing, and coherent imaging.

We examine the impact of top-of-screen promotions on viewing time at ABEMA, a leading video streaming platform in Japan. To this end, we conduct a large-scale randomized controlled trial. Given the non-standard distribution of user viewing times, we estimate distributional treatment effects. Our estimation results document that spotlighting content through these promotions effectively boosts user engagement across diverse content types. Notably, promoting short content proves most effective in that it not only retains users but also motivates them to watch subsequent episodes.

Credit risk models are a critical decision-support tool for financial institutions, yet tightening data-protection rules (e.g., GDPR, CCPA) increasingly prohibit cross-border sharing of borrower data, even as these models benefit from cross-institution learning. Traditional default prediction suffers from two limitations: binary classification ignores default timing, treating early defaulters (high loss) equivalently to late defaulters (low loss), and centralized training violates emerging regulatory constraints. We propose a Federated Survival Learning framework with Bayesian Differential Privacy (FSL-BDP) that models time-to-default trajectories without centralizing sensitive data. The framework provides Bayesian (data-dependent) differential privacy (DP) guarantees while enabling institutions to jointly learn risk dynamics. Experiments on three real-world credit datasets (LendingClub, SBA, Bondora) show that federation fundamentally alters the relative effectiveness of privacy mechanisms. While classical DP performs better than Bayesian DP in centralized settings, the latter benefits substantially more from federation (+7.0% vs +1.4%), achieving near parity of non-private performance and outperforming classical DP in the majority of participating clients. This ranking reversal yields a key decision-support insight: privacy mechanism selection should be evaluated in the target deployment architecture, rather than centralized benchmarks. These findings provide actionable guidance for practitioners designing privacy-preserving decision support systems in regulated, multi-institutional environments.

Maronna's and Tyler's $M$-estimators are among the most widely used robust estimators for scatter matrices. However, when the dimension of observations is relatively high, their performance can substantially deteriorate in certain situations, particularly in the presence of clustered outliers. To address this issue, we propose an estimator that shrinks the estimated precision matrix toward the identity matrix. We derive a sufficient condition for its existence, discuss its statistical interpretation, and establish upper and lower bounds for its additive finite sample breakdown point. Numerical experiments confirm the robustness of the proposed method.

In this paper, we propose a distributed framework for reducing the dimensionality of high-dimensional, large-scale, heterogeneous matrix-variate time series data using a factor model. The data are first partitioned column-wise (or row-wise) and allocated to node servers, where each node estimates the row (or column) loading matrix via two-dimensional tensor PCA. These local estimates are then transmitted to a central server and aggregated, followed by a final PCA step to obtain the global row (or column) loading matrix estimator. Given the estimated loading matrices, the corresponding factor matrices are subsequently computed. Unlike existing distributed approaches, our framework preserves the latent matrix structure, thereby improving computational efficiency and enhancing information utilization. We also discuss row- and column-wise clustering procedures for settings in which the group memberships are unknown. Furthermore, we extend the analysis to unit-root nonstationary matrix-variate time series. Asymptotic properties of the proposed method are derived for the diverging dimension of the data in each computing unit and the sample size $T$. Simulation results assess the computational efficiency and estimation accuracy of the proposed framework, and real data applications further validate its predictive performance.

We introduce a Markov chain Monte Carlo algorithm based on Sub-Cauchy Projection, a geometric transformation that generalizes stereographic projection by mapping Euclidean space into a spherical cap of a hyper-sphere, referred to as the complement of the dark side of the moon. We prove that our proposed method is uniformly ergodic for sub-Cauchy targets, namely targets whose tails are at most as heavy as a multidimensional Cauchy distribution, and show empirically its performance for challenging high-dimensional problems. The simplicity and broad applicability of our approach open new opportunities for Bayesian modeling and computation with heavy-tailed distributions in settings where most existing methods are unreliable.

Classical approaches to penalized smoothing often rely on basis or kernel expansions, which constrain the estimator to a fixed span and may be restrictive for discretely observed data. We instead regularize a single noisy trajectory directly on its observation grid, using difference operators that remain genuine finite-difference approximations to derivatives while being statistically normalized and mutually decorrelated under a reference noise law. We extend this white-noise construction to a parametric family of covariance-adapted reference geometries, with the parameter estimated by generalized method of moments from an independent calibration sample. A first-order plug-in expansion shows that, within this moment family, efficient calibration minimizes the leading-order discrepancy between the plug-in and oracle covariance-adapted smoothers. Numerical experiments confirm the predicted rate at which this plug-in discrepancy vanishes and compare the resulting reconstruction against conventional discrete, basis, and kernel smoothers.

We propose a framework for contextual distributionally robust optimization (DRO) that considers the causal and continuous structure of the underlying distribution, and we develop an interpretable and tractable decision rule. We first introduce the causal Sinkhorn discrepancy (CSD), an entropy-regularized causal Wasserstein distance that encourages continuous transport plans while preserving causal consistency. We then formulate a contextual DRO model with a CSD-based ambiguity set, termed Causal Sinkhorn DRO (Causal-SDRO), and derive its strong dual reformulation, where the worst-case distribution is characterized as a mixture of Gibbs distributions. To obtain an (infinite-dimensional) optimal policy, we propose a soft regression forest (SRF) decision rule: it preserves the interpretability of classical decision trees while being fully parametric, differentiable, and Lipschitz-smooth, enabling intrinsic interpretation from both global and local perspectives. To solve the Causal-SDRO with parametric decision rules, we develop an efficient stochastic compositional gradient algorithm that converges to an $\varepsilon$-stationary point at a rate of $\mathcal{O}(\varepsilon^{-4})$, matching that of standard stochastic gradient descent. Finally, we validate our method through numerical experiments on synthetic and real-world datasets, demonstrating its superior performance and interpretability.

Physical activity is crucial for human health. With the increasing availability of large-scale mobile health data, strong associations have been found between physical activity and various diseases. However, accurately capturing this complex relationship is challenging, possibly because it varies across different subgroups of subjects, especially in large-scale datasets. To fill this gap, we propose a generalized heterogeneous functional method which simultaneously estimates functional effects and identifies subgroups within the generalized functional regression framework. The proposed method captures subgroup-specific functional relationships between physical activity and diseases, providing a more nuanced understanding of these associations. Additionally, we develop a pre-clustering method that enhances computational efficiency for large-scale data through a finer partition of subjects compared to true subgroups. We further introduce a testing procedure to assess whether the different subgroups exhibit distinct functional effects. In the real data application, we examine the impact of physical activity on the risk of dementia using the UK Biobank dataset, which includes over 96,433 participants. Our proposed method outperforms existing methods in future-day prediction accuracy, identifying three distinct subgroups, with detailed scientific interpretations for each subgroup. We also demonstrate the theoretical consistency of our methods. Codes implementing the proposed method are available at: https://github.com/xiaojing777/GHFM.

Outlier detection (OD) aims to identify abnormal instances, known as outliers or anomalies, by learning typical patterns of normal data, or inliers. Performing OD under an unsupervised regime-without any information about anomalous instances in the training data-is challenging. A recently observed phenomenon, known as the inlier-memorization (IM) effect, where deep generative models (DGMs) tend to memorize inlier patterns during early training, provides a promising signal for distinguishing outliers. However, existing unsupervised approaches that rely solely on the IM effect still struggle when inliers and outliers are not well-separated or when outliers form dense clusters. To address these limitations, we incorporate active learning to selectively acquire informative labels, and propose IMBoost, a novel framework that explicitly reinforces the IM effect to improve outlier detection. Our method consists of two stages: 1) a warm-up phase that induces and promotes the IM effect, and 2) a polarization phase in which actively queried samples are used to maximize the discrepancy between inlier and outlier scores. In particular, we propose a novel query strategy and tailored loss function in the polarization phase to effectively identify informative samples and fully leverage the limited labeling budget. We provide a theoretical analysis showing that the IMBoost consistently decreases inlier risk while increasing outlier risk throughout training, thereby amplifying their separation. Extensive experiments on diverse benchmark datasets demonstrate that IMBoost not only significantly outperforms state-of-the-art active OD methods but also requires substantially less computational cost.

Constant rescaling of a Riemannian metric appears in many computational settings, often through a global scale parameter that is introduced either explicitly or implicitly. Although this operation is elementary, its consequences are not always made clear in practice and may be confused with changes in curvature, manifold structure, or coordinate representation. In this note we provide a short, self-contained account of constant metric scaling on arbitrary Riemannian manifolds. We distinguish between quantities that change under such a scaling, including norms, distances, volume elements, and gradient magnitudes, and geometric objects that remain invariant, such as the Levi–Civita connection, geodesics, exponential and logarithmic maps, and parallel transport. We also discuss implications for Riemannian optimization, where constant metric scaling can often be interpreted as a global rescaling of step sizes rather than a modification of the underlying geometry. The goal of this note is purely expository and is intended to clarify how a global metric scale parameter can be introduced in Riemannian computation without altering the geometric structures on which these methods rely.

This study uses connected vehicle data to analyze speeding behavior on residential roads. A scalable pipeline processes trajectory data and supplements missing speed limits to generate summaries at OpenStreetMap's way ID level. The findings reveal a highly skewed distribution of both aggressive and reckless speeding. Based on a case study of Charlottesville, VA's connected vehicle data on residential roads, we found that 38% of segments had at least one instance of aggressive speeding, and 20% had at least one instance of reckless speeding. In addition, night time speeding is 27 times more prevalent than day time, and extreme violations on specific road segments highlight how severe the issue can be. Several segments rank among the top 10 for both aggressive and reckless speedings, indicating that there exist high-risk residential roads. These findings support the need for both spatial and behavioral interventions. The analysis provides a rich foundation for policy and planning, offering a valuable complement to traditional enforcement and planning tools. In conclusion, this framework sets the foundation for future applications in traffic safety analytics, demonstrating the growing potential of telematics data to inform safer, more livable communities.