Papers

We introduce BayesChange, a computationally efficient R package, built on C++, for Bayesian change point detection and clustering of observations sharing common change points. While many R packages exist for change point analysis, BayesChange offers methods not currently available elsewhere. The core functions are implemented in C++ to ensures computational efficiency, while an R user interface simplifies the package usage. The BayesChange package includes two R wrappers that integrate the C++ backend functions, along with S3 methods for summarizing the results. We present the theory beyond each method, the algorithms for posterior simulation and we illustrate the package's usage through synthetic examples.

Heavy-tailed phenomena appear across diverse domains –from wealth and firm sizes in economics to network traffic, biological systems, and physical processes– characterized by the disproportionate influence of extreme values. These distributions challenge classical statistical models, as their tails decay too slowly for conventional approximations to hold. Among their key descriptive measures are quantile contributions, which quantify the proportion of a total quantity (such as income, energy, or risk) attributed to observations above a given quantile threshold. This paper presents a theoretical study of the quantile contribution statistic and its relationship with order statistics. We derive a closed-form expression for the joint cumulative distribution function (CDF) of order statistics and, based on it, obtain an explicit CDF for quantile contributions applicable to small samples. We then investigate the asymptotic behavior of these contributions as the sample size increases, establishing the asymptotic normality of the numerator and characterizing the limiting distribution of the quantile contribution. Finally, simulation studies illustrate the convergence properties and empirical accuracy of the theoretical results, providing a foundation for applying quantile contributions in the analysis of heavy-tailed data.

A sum-wise formulation is proposed for the Kaplan-Meier product limit estimator of partially right-censored survival data. The population estimator is expressed as a sum over the individual units' empirical and semi-empirical contributions, for observed and censored failures, respectively. This intuitive decomposition is applied to visualize the different contributions of failed and censored units to the overall population estimator.

A fundamental challenge in developing general learning algorithms is their tendency to forget past knowledge as they adapt to new data. Addressing this problem requires a principled understanding of forgetting. Yet, despite decades of study, no unified definition has emerged that offers insight into the underlying dynamics of learning. We propose an algorithm- and task-agnostic theory that characterises forgetting as a lack of self-consistency in a learner's predictive distribution, manifesting as a loss of predictive information. Our theory naturally yields a general measure of an algorithm's propensity to forget, proves that exact Bayesian inference allows for adaptation without forgetting, and provides a tautological explanation for why generative models forget when trained on their own synthetic outputs. To validate these claims, we design a comprehensive set of experiments that span classification, regression, generative modelling, and reinforcement learning. We demonstrate that forgetting is present across all deep learning settings and plays a significant role in determining learning efficiency. Together, these results establish a principled understanding of forgetting and lay the foundation for analysing and improving the information retention capabilities of general learning algorithms.

Multi-site studies are widely used across the social sciences, education, and medicine to learn the effects of interventions across heterogeneous populations. But deployment populations, on which a given intervention is implemented, may differ systematically from the experimental population on which the intervention is evaluated. I propose a method for selecting sites with minimax guarantees against directional distribution shifts, which allows the researcher to guard against specified threats to the external validity of a multi-site study. I use Wasserstein Distributionally Robust Optimization to select sites with guarantees against directional shifts, provide a novel cutting-plane algorithm to implement the approach, and prove its optimality with respect to minimizing worst-case shifts. I highlight the empirical performance of the method with with two semi-synthetic simulations: I re-analyse the Tennessee STAR experiment, under the assumption that the study population underrepresented low-SES schools, and \citet{Crepon_2015}, which assesses the impact of increasing access to microcredit in rural Morocco, under the assumption that the study population was skewed towards villages with larger market potential. In these examples, I show that the WDRO method minimizes mean bias, average mean squared error, and $W_1$ distance under specified directional distribution shifts.

The application of causal discovery to diseases like Alzheimer's (AD) is limited by the static graph assumptions of most methods; such models cannot account for an evolving pathophysiology, modulated by a latent disease pseudotime. We propose to apply an existing latent variable model to real-world AD data, inferring a pseudotime that orders patients along a data-driven disease trajectory independent of chronological age, then learning how causal relationships evolve. Pseudotime outperformed age in predicting diagnosis (AUC 0.82 vs 0.59). Incorporating minimal, disease-agnostic background knowledge substantially improved graph accuracy and orientation. Our framework reveals dynamic interactions between novel (NfL, GFAP) and established AD markers, enabling practical causal discovery despite violated assumptions.

Resilient supply chains are critical, especially for Original Equipment Manufacturers (OEMs) that power today's digital economy. Safety Stock dimensioning-the computation of the appropriate safety stock quantity-is one of several mechanisms to ensure supply chain resiliency, as it protects the supply chain against demand and supply uncertainties. Unfortunately, the major approaches to dimensioning safety stock heavily assume that demand is normally distributed and ignore future demand variability, limiting their applicability in manufacturing contexts where demand is non-normal, intermittent, and highly skewed. In this paper, we propose a data-driven approach that relaxes the assumption of normality, enabling the demand distribution of each inventory item to be analytically determined using Kernel Density Estimation. Also, we extended the analysis from historical demand variability to forecasted demand variability. We evaluated the proposed approach against a normal distribution model in a near-world inventory replenishment simulation. Afterwards, we used a linear optimization model to determine the optimal safety stock configuration. The results from the simulation and linear optimization models showed that the data-driven approach outperformed traditional approaches. In particular, the data-driven approach achieved the desired service levels at lower safety stock levels than the conventional approaches.

M-estimation, aka empirical risk minimization, is at the heart of statistics and machine learning: Classification, regression, location estimation, etc. Asymptotic theory is well understood when the loss satisfies some smoothness assumptions and its derivatives are dominated locally. However, these conditions are typically technical and can be too restrictive or heavy to check. Here, we consider the case of a convex loss function, which may not even be differentiable: We establish an asymptotic theory for M-estimation with convex loss (which needs not be differentiable) under convex constraints. We show that the asymptotic distributions of the corresponding M-estimators depend on an interplay between the loss function and the boundary structure of the set of constraints. We extend our results to U-estimators, building on the asymptotic theory of U-statistics. Applications of our work include, among other, robust location/scatter estimation, estimation of deepest points relative to depth functions such as Oja's depth, etc.

Understanding feature-outcome associations in high-dimensional data remains challenging when relationships vary across subpopulations, yet standard methods assuming global associations miss context-dependent patterns, reducing statistical power and interpretability. We develop a geometric decomposition framework offering two strategies for partitioning inference problems into regional analyses on data-derived Riemannian graphs. Gradient flow decomposition uses path-monotonicity-validated discrete Morse theory to partition samples into gradient flow cells where outcomes exhibit monotonic behavior. Co-monotonicity decomposition utilizes vertex-level coefficients that provide context-dependent versions of the classical Pearson correlation: these coefficients measure edge-based directional concordance between outcome and features, or between feature pairs, defining embeddings of samples into association space. These embeddings induce Riemannian k-NN graphs on which biclustering identifies co-monotonicity cells (coherent regions) and feature modules. This extends naturally to multi-modal integration across multiple feature sets. Both strategies apply independently or jointly, with Bayesian posterior sampling providing credible intervals.

In the theory of optimal transport, the Knothe-Rosenblatt (KR) rearrangement provides an explicit construction to map between two probability measures by building one-dimensional transformations from the marginal conditionals of one measure to the other. The KR map has shown to be useful in different realms of mathematics and statistics, from proving functional inequalities to designing methodologies for sampling conditional distributions. It is known that the KR rearrangement can be obtained as the limit of a sequence of optimal transport maps with a weighted quadratic cost. We extend these results in this work by showing that one can obtain the KR map as a limit of maps that solve a relaxation of the weighted-cost optimal transport problem with a soft-constraint for the target distribution. In addition, we show that this procedure also applies to the construction of triangular velocity fields via dynamic optimal transport yielding optimal velocity fields. This justifies various variational methodologies for estimating KR maps in practice by minimizing a divergence between the target and pushforward measure through an approximate map. Moreover, it opens the possibilities for novel static and dynamic OT estimators for KR maps.

PDEs arise ubiquitously in science and engineering, where solutions depend on parameters (physical properties, boundary conditions, geometry). Traditional numerical methods require re-solving the PDE for each parameter, making parameter space exploration prohibitively expensive. Recent machine learning advances, particularly physics-informed neural networks (PINNs) and neural operators, have revolutionized parametric PDE solving by learning solution operators that generalize across parameter spaces. We critically analyze two main paradigms: (1) PINNs, which embed physical laws as soft constraints and excel at inverse problems with sparse data, and (2) neural operators (e.g., DeepONet, Fourier Neural Operator), which learn mappings between infinite-dimensional function spaces and achieve unprecedented generalization. Through comparisons across fluid dynamics, solid mechanics, heat transfer, and electromagnetics, we show neural operators can achieve computational speedups of $10^3$ to $10^5$ times faster than traditional solvers for multi-query scenarios, while maintaining comparable accuracy. We provide practical guidance for method selection, discuss theoretical foundations (universal approximation, convergence), and identify critical open challenges: high-dimensional parameters, complex geometries, and out-of-distribution generalization. This work establishes a unified framework for understanding parametric PDE solvers via operator learning, offering a comprehensive, incrementally updated resource for this rapidly evolving field

This study clarifies the relationship between Riesz regression [Chernozhukov et al., 2021] and density ratio estimation (DRE) in causal inference problems, such as average treatment effect estimation. We first show that the Riesz representer can be written as a signed density ratio and then demonstrate that the Riesz regression objective coincides with the least-squares importance fitting criterion [Kanamori et al., 2009]. Although Riesz regression applies to a broad class of representer estimation problems, this equivalence with DRE allows us to transfer existing DRE results, including convergence rate analyses, generalizations based on Bregman divergence minimization, and regularization techniques for flexible models such as neural networks.

In this paper, we analyze the asymptotic behavior of a system of interacting reinforced stochastic processes $({\bf Z}_n, {\bf N}_n)_n$ on a directed network of $N$ agents. The system is defined by the coupled dynamics ${\bf Z}_{n+1}=(1-r_{n}){\bf Z}_{n}+r_{n}{\bf X}_{n+1}$ and ${\bf N}_{n+1}=(1-\frac{1}{n+1}){\bf N}_n+\frac{1}{n+1}{\bf X}_{n+1}$, where agent actions $\mathbb{P}(X_{n+1,j}=1\mid{\cal F}_n)=\sum_{h} w_{hj}Z_{nh}$ are governed by a column-normalized adjacency matrix ${\bf W}$, and $r_n \sim cn^{-\gamma}$ with $\gamma \in (1/2, 1]$. Existing asymptotic theory has largely been restricted to irreducible and diagonalizable ${\bf W}$. We extend this analysis to the broader and more practical class of reducible and non-diagonalizable matrices ${\bf W}$ possessing a block upper-triangular form, which models hierarchical influence. We first establish synchronization, proving $({\bf Z}^\top_n, {\bf N}^\top_n)^\top \to Z_\infty {\bf 1}$ almost surely, where the distribution of the limit $Z_\infty$ is shown to be determined solely by the internal dynamics of the leading subgroup. Furthermore, we establish a joint central limit theorem for $({\bf Z}_n,{\bf N}_n)_n$, revealing how the spectral properties and Jordan block structure of ${\bf W}$ govern second-order fluctuations. We demonstrate that the convergence rates and the limiting covariance structure exhibit a phase transition dependent on $\gamma$ and the spectral properties of ${\bf W}$. Crucially, we explicitly characterize how the non-diagonalizability of ${\bf W}$ fundamentally alters the asymptotic covariance and introduces new logarithmic scaling factors in the critical case ($\gamma=1$). These results provide a probabilistic foundation for statistical inference on such hierarchical network structures.

Generative Bayesian Filtering (GBF) provides a powerful and flexible framework for performing posterior inference in complex nonlinear and non-Gaussian state-space models. Our approach extends Generative Bayesian Computation (GBC) to dynamic settings, enabling recursive posterior inference using simulation-based methods powered by deep neural networks. GBF does not require explicit density evaluations, making it particularly effective when observation or transition distributions are analytically intractable. To address parameter learning, we introduce the Generative-Gibbs sampler, which bypasses explicit density evaluation by iteratively sampling each variable from its implicit full conditional distribution. Such technique is broadly applicable and enables inference in hierarchical Bayesian models with intractable densities, including state-space models. We assess the performance of the proposed methodologies through both simulated and empirical studies, including the estimation of $\alpha$-stable stochastic volatility models. Our findings indicate that GBF significantly outperforms existing likelihood-free approaches in accuracy and robustness when dealing with intractable state-space models.

We present a new benchmarking study comparing a boundary-constrained Ehrenpreis–Palamodov Gaussian Process (B-EPGP) surrogate with a classical finite element method combined with Crank–Nicolson time stepping (CN-FEM) for solving the two-dimensional wave equation with homogeneous Dirichlet boundary conditions. The B-EPGP construction leverages exponential-polynomial bases derived from the characteristic variety to enforce the PDE and boundary conditions exactly and employs penalized least squares to estimate the coefficients. To ensure fairness across paradigms, we introduce a degrees-of-freedom (DoF) matching protocol. Under matched DoF, B-EPGP consistently attains lower space-time $L^2$-error and maximum-in-time $L^{2}$-error in space than CN-FEM, improving accuracy by roughly two orders of magnitude.

We propose a unified class of calibration weighting methods based on weighted generalized entropy to handle missing at random (MAR) data with improved stability and efficiency. The proposed generalized entropy calibration (GEC) formulates weight construction as a convex optimization program that unifies entropy-based approaches and generalized regression weighting. Double robustness is achieved by augmenting standard covariate balancing with a debiasing constraint tied to the propensity score model and a Neyman-orthogonal constraint that removes first-order sensitivity to nuisance estimation. Selection of the weights on the entropy function can lead to the optimal calibration estimator under a correctly specified outcome regression model. The proposed GEC weighting ha a nice geometric characterization: the GEC solution is the Bregman projection of the initial weights onto a constraint set, which yields a generalized Pythagorean identity and a nested decomposition that quantifies the incremental distance paid for additional constraints. We also develop a high-dimensional extension with soft calibration and a projection calibration constraint that preserves doubly robust inference. Two simulation studies are presented to compare the performance of the proposed method with the existing methods.

Market generators using deep generative models have shown promise for synthetic financial data generation, but existing approaches lack causal reasoning capabilities essential for counterfactual analysis and risk assessment. We propose a Time-series Neural Causal Model VAE (TNCM-VAE) that combines variational autoencoders with structural causal models to generate counterfactual financial time series while preserving both temporal dependencies and causal relationships. Our approach enforces causal constraints through directed acyclic graphs in the decoder architecture and employs the causal Wasserstein distance for training. We validate our method on synthetic autoregressive models inspired by the Ornstein-Uhlenbeck process, demonstrating superior performance in counterfactual probability estimation with L1 distances as low as 0.03-0.10 compared to ground truth. The model enables financial stress testing, scenario analysis, and enhanced backtesting by generating plausible counterfactual market trajectories that respect underlying causal mechanisms.

We consider the problem of testing for differences in group-specific slopes between the selected groups in panel data identified via k-means clustering. In this setting, the classical Wald-type test statistic is problematic because it produces an extremely inflated type I error probability. The underlying reason is that the same dataset is used to identify the group structure and construct the test statistic, simultaneously. This creates dependence between the selection and inference stages. To address this issue, we propose a valid selective inference approach conditional on the selection event to account for the selection effect. We formally define the selective type I error and describe how to efficiently compute the correct p-values for clusters obtained using k-means clustering. Furthermore, the same idea can be extended to test for differences in coefficients due to a single covariate and can be incorporated into the GMM estimation framework. Simulation studies show that our method has satisfactory finite sample performance. We apply this method to explore the heterogeneous relationships between economic growth and the $CO_2$ emission across countries for which some new findings are discovered. An R package TestHomoPanel is provided to implement the proposed selective inference framework for panel data.

Positron emission tomography (PET) imaging is widely used in a number of clinical applications, including cancer and Alzheimer's disease (AD) diagnosis, monitoring of disease development, and treatment effect evaluation. Statistical modeling of PET imaging is essential to address continually emerging scientific questions in these research fields, including hypotheses related to evaluation of effects of disease modifying treatments on amyloid reduction in AD and associations between amyloid reduction and cognitive function, among many others. In this paper, we provide background information and tools for statisticians interested in developing statistical models for PET imaging to pre-process and prepare data for analysis. We introduce our novel pre-processing and visualization tool TRAECR (Template registration, MRI-PET co-Registration, Anatomical brain Extraction and COMBAT/RAVEL harmonization) to facilitate data preparation for statistical analysis.

We study ranking and selection under input uncertainty in settings where additional data cannot be collected. We propose the Nonparametric Input-Output Uncertainty Comparisons (NIOU-C) procedure to construct a confidence set that includes the optimal solution with a user-specified probability. We construct an ambiguity set of input distributions using empirical likelihood and approximate the mean performance of each solution using a linear functional representation of the input distributions. By solving optimization problems evaluating worst-case pairwise mean differences within the ambiguity set, we build a confidence set of solutions indistinguishable from the optimum. We characterize sample size requirements for NIOU-C to achieve the asymptotic validity under mild conditions. Moreover, we propose an extension to NIOU-C, NIOU-C:E, that mitigates conservatism and yields a smaller confidence set. In numerical experiments, NIOU-C provides a smaller confidence set that includes the optimum more frequently than a parametric procedure that takes advantage of the parametric distribution families.

Time series forecasting is essential across domains from finance to supply chain management. This paper introduces ForecastGAN, a novel decomposition based adversarial framework addressing limitations in existing approaches for multi-horizon predictions. Although transformer models excel in long-term forecasting, they often underperform in short-term scenarios and typically ignore categorical features. ForecastGAN operates through three integrated modules: a Decomposition Module that extracts seasonality and trend components; a Model Selection Module that identifies optimal neural network configurations based on forecasting horizon; and an Adversarial Training Module that enhances prediction robustness through Conditional Generative Adversarial Network training. Unlike conventional approaches, ForecastGAN effectively integrates both numerical and categorical features. We validate our framework on eleven benchmark multivariate time series datasets that span various forecasting horizons. The results show that ForecastGAN consistently outperforms state-of-the-art transformer models for short-term forecasting while remaining competitive for long-term horizons. This research establishes a more generalizable approach to time series forecasting that adapts to specific contexts while maintaining strong performance across diverse data characteristics without extensive hyperparameter tuning.

The design of informatively rich input signals is essential for accurate system identification, yet classical Fisher-information-based methods are inherently local and often inadequate in the presence of significant model uncertainty and nonlinearity. This paper develops a Bayesian approach that uses the mutual information (MI) between observations and parameters as the utility function. To address the computational intractability of the MI, we maximize a tractable MI lower bound. The method is then applied to the design of an input signals for the identification of quasi-linear stochastic dynamical systems. Evaluating the MI lower bound requires inversion of large covariance matrices whose dimensions scale with the number of data points $N$. To overcome this problem, an algorithm that reduces the dimension of the matrices to be inverted by a factor of $N$ is developed, making the approach feasible for long experiments. The proposed Bayesian method is compared with the average D-optimal design method, a semi-Bayesian approach, and its advantages are demonstrated. The effectiveness of the proposed method is further illustrated through four examples, including atomic sensor models, where the input signals that generates large MI are especially important for reducing the estimation error.

Bayesian experimental design (BED) provides a principled framework for optimizing data collection by choosing experiments that are maximally informative about unknown parameters. However, existing methods cannot deal with the joint challenge of (a) partially observable dynamical systems, where only noisy and incomplete observations are available, and (b) fully online inference, which updates posterior distributions and selects designs sequentially in a computationally efficient manner. Under partial observability, dynamical systems are naturally modeled as state-space models (SSMs), where latent states mediate the link between parameters and data, making the likelihood – and thus information-theoretic objectives like the expected information gain (EIG) – intractable. We address these challenges by deriving new estimators of the EIG and its gradient that explicitly marginalize latent states, enabling scalable stochastic optimization in nonlinear SSMs. Our approach leverages nested particle filters for efficient online state-parameter inference with convergence guarantees. Applications to realistic models, such as the susceptible-infectious-recovered (SIR) and a moving source location task, show that our framework successfully handles both partial observability and online inference.

Energy markets exhibit complex causal relationships between weather patterns, generation technologies, and price formation, with regime changes occurring continuously rather than at discrete break points. Current approaches model electricity prices without explicit causal interpretation or counterfactual reasoning capabilities. We introduce Augmented Time Series Causal Models (ATSCM) for energy markets, extending counterfactual reasoning frameworks to multivariate temporal data with learned causal structure. Our approach models energy systems through interpretable factors (weather, generation mix, demand patterns), rich grid dynamics, and observable market variables. We integrate neural causal discovery to learn time-varying causal graphs without requiring ground truth DAGs. Applied to real-world electricity price data, ATSCM enables novel counterfactual queries such as "What would prices be under different renewable generation scenarios?".

We establish strong connections between two fundamental nonlinear 0/1 optimization problems coming from the area of experimental design, namely maximum-entropy sampling and 0/1 D-optimality. The connections are based on maps between instances, and we analyze the behavior of these maps. Using these maps, we transport basic upper-bounding methods between these two problems. Further, we establish results relating how different branch-and-bound schemes based on these maps compare. Additionally, we observe some surprising numerical results, where bounding methods that did not seem promising in their direct application to real-data MESP instances, are now useful for MESP instances that come from 0/1 D-optimality.

This paper introduces a matrix-variate regression model for analyzing multivariate data observed across spatial locations and over time. The model's design incorporates a mean structure that links covariates to the response matrix and a separable covariance structure, based on a Kronecker product, to capture spatial and temporal dependencies efficiently. We derive maximum likelihood estimators for all model parameters. A simulation study validates the model, showing its effectiveness in parameter recovery across different spatial resolutions. Finally, an application to real-world data on agricultural and livestock production from Brazilian municipalities showcases the model's practical utility in revealing structured spatio-temporal patterns of variation and covariate effects.

A central part of geometric statistics is to compute the Fréchet mean. This is a well-known intrinsic mean on a Riemannian manifold that minimizes the sum of squared Riemannian distances from the mean point to all other data points. The Fréchet mean is simple to define and generalizes the Euclidean mean, but for most manifolds even minimizing the Riemannian distance involves solving an optimization problem. Therefore, numerical computations of the Fréchet mean require solving an embedded optimization problem in each iteration. We introduce the GEORCE-FM algorithm to simultaneously compute the Fréchet mean and Riemannian distances in each iteration in a local chart, making it faster than previous methods. We extend the algorithm to Finsler manifolds and introduce an adaptive extension such that GEORCE-FM scales to a large number of data points. Theoretically, we show that GEORCE-FM has global convergence and local quadratic convergence and prove that the adaptive extension converges in expectation to the Fréchet mean. We further empirically demonstrate that GEORCE-FM outperforms existing baseline methods to estimate the Fréchet mean in terms of both accuracy and runtime.

This paper presents different recursive formulas for computing the marginals and the normalizing constant of a Gibbs distribution $\pi$: The common thread is the use of the underlying Markov properties of such processes. The procedures are illustrated with several examples, particularly the Ising model.

Minimum-volume nonnegative matrix factorization (min-vol NMF) has been used successfully in many applications, such as hyperspectral imaging, chemical kinetics, spectroscopy, topic modeling, and audio source separation. However, its robustness to noise has been a long-standing open problem. In this paper, we prove that min-vol NMF identifies the groundtruth factors in the presence of noise under a condition referred to as the expanded sufficiently scattered condition which requires the data points to be sufficiently well scattered in the latent simplex generated by the basis vectors.

Conformal prediction has emerged as a powerful framework for constructing distribution-free prediction sets with guaranteed coverage assuming only the exchangeability assumption. However, this assumption is often violated in online environments where data distributions evolve over time. Several recent approaches have been proposed to address this limitation, but, typically, they slowly adapt to distribution shifts because they update predictions only in a forward manner, that is, they generate a prediction for a newly observed data point while previously computed predictions are not updated. In this paper, we propose a novel online conformal inference method with retrospective adjustment, which is designed to achieve faster adaptation to distributional shifts. Our method leverages regression approaches with efficient leave-one-out update formulas to retroactively adjust past predictions when new data arrive, thereby aligning the entire set of predictions with the most recent data distribution. Through extensive numerical studies performed on both synthetic and real-world data sets, we show that the proposed approach achieves coverage close to the nominal level while reducing predictive interval width by up to approximately 30% compared to existing online conformal prediction methods, demonstrating improved statistical efficiency alongside faster adaptation.