2026-10-14 | | Total: 218
Simulation-based Bayesian inference for Chemical Reaction Networks (CRNs) must simultaneously solve two computationally demanding problems: discriminating between overlapping reaction pathways and estimating reaction rate constants. Classical approaches, such as Sequential Monte Carlo ABC, rely on handcrafted summary statistics and scale poorly across large parameter spaces. To overcome this, we introduce a neural inference architecture. Our pipeline utilizes shared 1D-CNNs and a permutation-invariant Set Attention mechanism to learn low-dimensional summary statistics from variable-length, multi-experiment time-series data. These statistics simultaneously condition a model selection network to compute posterior model probabilities and model-specific conditional Invertible Neural Networks (cINNs) that estimate parameter posterior distributions. We demonstrate the practical utility of our approach on a real protodeboronation example, illustrating how it handles structural ambiguity and extracts kinetic parameters from competing catalyst deactivation pathways.
Mortgage credit risk may exhibit dependence across economically connected areas, even when they are geographically distant. However, current models typically represent cross-area dependence through fixed geographical neighbourhoods, potentially overlooking economic connections that extend beyond physical proximity and evolve over time. These models may fail to capture shared risks among geographically distant but economically similar areas, and cannot describe a changing dependence of risk; it therefore risks retaining economically obsolete connections while omitting newly emerging ones. This paper uses similarities of industrial composition as a measure of economic proximity and develops a discrete-time survival model with dynamic industrial graphs to capture evolving cross-area dependence in mortgage default risk. The resulting graphs are reconstructed annually, allowing both their topology and edge weights to change as regional industrial compositions evolve. Jointly combining a proper conditional autoregressive process with a first-order autoregressive process governs the temporal cross-area dependence and the temporal persistence in latent risk simultaneously. The model is estimated using Integrated Nested Laplace Approximation and evaluated on U.S. mortgage loan-level data against the static industrial graph specifications. The results show that the proposed dynamic model achieves better in-sample fit and out-of-sample discrimination, while its estimated parameters indicate substantial graph-induced cross-area dependence and temporal persistence. These findings demonstrate the value of allowing economically defined neighbourhoods to evolve when modelling mortgage default risk.
Design calculations let researchers gauge, before data is collected, whether a planned study is likely to detect an effect. In this paper, we develop such calculations for predictive models using ordinary least squares and ridge regression under leave-one-out cross-validation. To do so, we derive closed-form expressions relating the expected out-of-sample $R^2$ and the expected correlation between observed and predicted outcomes to four design-stage quantities: sample size, number of predictors, an assumed population $R^2$ ($ρ^2_{in}$) and a design efficiency factor $D^*$ (which summarizes how leverage is distributed across observations). Unlike the classical shrinkage formulas, which depend only on sample size and dimension, these are sensitive to how covariates are arranged -- to imbalance and skew. A simulation study assesses the accuracy of the formulas against realized values under a variety of data generating mechanisms. We apply the methods to a published study predicting personality from smartphone behavior and show that the design formulas would have indicated, before any analyses were run, that least squares could not be expected to beat the sample mean while ridge could. An R package, predpower, implements the calculations.
Impact-echo testing is a nondestructive testing method used to detect internal defects in concrete structures. However, boundary reflections, structured interference, and measurement noise can produce complex Fourier amplitude spectra, making defect classification difficult. We constructed a two-dimensional frequency--amplitude persistence image (FA-PI) from the Fourier amplitude spectrum. Using a fixed parameterization of exact multi-parameter persistent homology, each retained Fourier component is represented with frequency-bin location as one coordinate and amplitude as the other, and Gaussian kernels are then applied in this two-dimensional space. This construction has a direct analytical connection to the corresponding one-dimensional Gaussian-smoothed Fourier amplitude spectrum. FA-PI was evaluated on an experimental backfill-void dataset and a public bridge-deck impact-echo dataset, and compared with raw and Gaussian-smoothed Fourier amplitude spectra using a linear classifier with a Softmax output and identical signal-level data splits. The evaluation included original signals and signals corrupted by Gaussian noise, sinusoidal interference, or mixed noise, with matching noise types and signal-to-noise ratios during training, validation, and testing. Performance was assessed using accuracy and macro-averaged precision, recall, and F1-score. Across the two datasets, FA-PI consistently outperformed the raw and Gaussian-smoothed Fourier baselines in the noisy-condition summaries. These findings indicate that representing Fourier spectral components in a two-dimensional frequency--amplitude space can improve impact-echo defect classification beyond frequency-axis smoothing alone.
We consider the problem of estimating the expectation of a function (observable) with respect to a high-dimensional probability distribution, under the assumption that a low-dimensional collective variable map captures all the directions of slow mixing for the distribution. More precisely, under the assumption that the fibers of the collective variable map satisfy a log-Sobolev inequality, we establish non-asymptotic relative-error guarantees for this problem using the eigenvector method for umbrella sampling (EMUS). Our analysis applies to both log-Lipschitz observables and indicator of convex sets, and gives complexity bounds that are polynomial in regularity parameters and inverse relative accuracy, and exponential in the dimension of the collective variable space. As a relative error guarantee, our bounds apply to rare event estimation. We consider Gaussian biasing windows defined on a low-dimensional collective variable space and analyze a modified sparse matrix version of EMUS. The bounds quantify the effect of the map through mixing of the conditional distribution and fluctuation of the observable along the fibers.
We construct confidence sets for the entire Fréchet mean set by inverting studentized pairwise risk inequalities. Gaussian multiplier calibration localized around empirical near minimizers gives asymptotic coverage at least the nominal level, uniformly over distribution classes satisfying empirical-process conditions. The construction requires neither unique ness nor a manifold structure. We characterize the limiting binding directions and give conditions for exact asymptotic coverage. On the intrinsic circle, exactness holds for every fixed non-point-mass law and uniformly over a family in which two means merge at a quar tic minimum. With ideal multiplier quantiles, the construction reproduces the Euclidean Wald set under both nonsingular and singular covariance. Under regular local conditions, it is first-order equivalent to intrinsic Wald inference.
We study the role of design smoothness in the minimax estimation of conditional average treatment effects at a fixed covariate point. For design densities bounded above and away from zero, without further regularity, we derive upper and lower risk bounds with an arbitrarily small gap in the convergence-rate exponent. Under uniform design, we obtain matching upper and lower bounds at interior points in the regime of low nuisance smoothness. These results show that rough covariate designs can lead to slower rates than uniform designs, even under the same smoothness assumptions on the treatment effect and nuisance functions, highlighting the role of the covariate distribution in pointwise causal estimation.
Every October, the Nobel Prizes in Physics, Chemistry, Physiology or Medicine and Economic Sciences are announced, with one to three prizes in each discipline. Of the 769 Laureates (since 1901) in these four disciplines only 4% are women (30 Laureates). In 2026, there are eight Laureates within these disciplines and they are all men. Based on historical data across four fields and a Bayesian hierarchical model, we update the 8 year old model on bias in Nobel Prizes and confirm the gender gap, with exceedingly large confidence.
CMinHash promises a new direction for designing MinHash algorithms with reduced variance and storage space via circulant permutations. However, the original variance analysis is difficult to extend to other MinHash estimators. We present a new approach that gives explicit and easily comparable approximate variance expressions, and use it to analyze the Minner estimator. Numerical experiments show that our variance expression matches the empirical mean square error (MSE) of estimates. Our work thereby demonstrates how to apply circulant permutation to other MinHash estimators and compute their approximate variance.
In many applications one object is chosen to represent a sample and the least typical objects are set aside, usually without any measure of uncertainty. We measure how central an object is by its reference depth, the probability that a random reference point is closer to it than to a random opponent. Using only pairwise dissimilarities, we construct simultaneous intervals for all ranks, a centre set that contains the most central observation under stated conditions, and outlier sets that contain only objects among the least central ones, each with asymptotic probability at least a prescribed level. These guarantees hold even though the ranked objects are themselves random. For symmetric distributions on the real line, the centre set has the smallest possible order of size up to a logarithmic factor, and in dimension three or more it reduces to the most central observation with probability tending to one. For typical meteorological years the centre set usually contains several candidate years, so the data do not determine the year on which building-energy simulations rest; for three collections of gene trees it contains one tree, whose features can be reported as those of the most central tree.
Sampling from high-dimensional, heavy-tailed distributions is a fundamental challenge in computational statistics, as many standard Markov chain Monte Carlo (MCMC) methods mix poorly in such settings. Recently, Stereographic MCMC [Yang et al., 2024] and the Sub-Cauchy Projection Sampler [Grazzi et al., 2026] have been shown to perform successfully on such tasks. However, establishing their non-asymptotic convergence properties remains an important open problem. In this work, we fill this gap by establishing non-asymptotic upper bounds on the mixing time of the stereographic projection and sub-Cauchy projection samplers. Our results demonstrate that, under certain conditions on the target and initial distributions, the mixing time is polynomial in dimension for a broad class of distributions, including light- and heavy-tailed cases. Motivated by the theoretical analysis, we further establish a new weighted isoperimetric inequality that extends the classical version for log-concave distributions to the heavy-tailed setting with optimal dimension dependence. The proof techniques provide new insights into the geometric properties of heavy-tailed distributions that govern rapid mixing in high dimensions.
Recently, Morris et al. introduced TinyLoRA, a strategy for LLM fine-tuning which introduces arbitrarily few parameters, and demonstrated accuracy improvements using less than $10$ parameters rather than the ${\gtrsim 10,000}$ typically demanded by LoRA-based methods. However, that study continued to use gradient-based local search methods of the kind deployed for training multi-billion parameter LLMs in the first place. This article is predicated on the observation that, in the low dimensional setting, global optimization methods become tractable, unlocking Large language Model Bayesian Optimization (LlMBO). LlMBO is sample-efficient, gradient-free, memory-light, parallelizes across devices without synchronization, and, unlike policy-gradient methods, can optimize greedy-decoding accuracy directly. We describe numerical experiments comparing LlMBO against GRPO on fine-tuning tasks for a variety of models on a math benchmark suite which show that LlMBO can be a suitable alternative in this regime, matching or exceeding GRPO's fine tuning gains, while using a fraction of the compute and memory.
In the past two decades, the Bayesian approach to inverse problems in infinite dimensional spaces has seen significant developments in theory, Markov chain Monte Carlo (MCMC) methods, and in applications. The Bayesian approach in finite dimension produces posterior distributions of finite dimensional random variables. In infinite dimensions, the desired objects are random fields and their laws as measures on Hilbert or Banach spaces. There are many textbooks and tutorials available that explain Bayesian methods and MCMC in finite dimensions. To enter the field of Bayesian inverse problems in infinite dimensions, a newcomer has to engage with a quite diverse and extensive literature. The purpose of this article is to provide a basic path to the subject that is easily accessible to a reader with general knowledge in probability and mathematical analysis. We start with a swift introduction to the finite dimensional case. The infinite dimensional case is introduced for Gaussian priors on Hilbert spaces and Gaussian-type likelihood functions. In this situation, mere continuity of the forward map suffices to prove that the posterior measure is well-defined. We also give some new and simple stability results. Two application examples are elaborated, one finite dimensional, one infinite dimensional.
Expected improvement (EI) based on stochastic kriging (SK) surrogate is widely used in continuous simulation optimization, but unknown heteroscedastic noise on simulation output makes three basic design choices consequential: the incumbent used to define improvement, the replication level at each sampled point, and the terminal recommendation. Existing methods typically either assume the noise levels are known or choose these components separately through heuristic rules or empirical calibration. We instead derive a high-probability upper bound on the simple regret and systematically choose the three components within general rule classes to tighten this bound. We first establish a finite-sample confidence bound for the SK predictor with unknown and estimated noise levels as a foundation for the three design choices. We then analyze general posterior-quantile families for both the incumbent and the terminal recommendation. For the incumbent, minimizing the regret bound selects a lower quantile coefficient that matches the coefficient in the SK predictor confidence bound. For the terminal recommendation, minimizing the same bound instead selects the sampled point with the best upper posterior quantile under the same coefficient. These opposite quantile choices reflect the different roles of search and implementation in simulation optimization. Intermediate simulations are exploratory with an optimistic incumbent, while the terminal recommendation is implemented in the real system and is therefore relatively conservative about the remaining posterior uncertainty. Under mild conditions satisfied by common kernels, minimizing a tractable kernel-specific bound yields $Θ(\log T)$ replications. To our knowledge, this is the first simple-regret guarantee for a continuous-domain EI-type simulation optimization method that learns unknown heteroscedastic variances through replications.
We propose a multilevel main effects matrix factor model (MMEFM) to extract meaningful row and column effects from groups of matrix data. As a generalization of a multilevel matrix factor model, we provide rigorous definitions of MMEFM and identifications, together with iterative algorithms for the estimators of all global and local components in the main effects and common component, including global/local main effects and core ranks, with rates of convergence spelt out. An important feature of the model is that all main effects can be nonstationary and still be consistently estimated at each time point. At the same time, our estimated global common factors are alignment independent, and hence are truly global, and can be presented and given interpretations more easily. We propose a Bootstrap procedure for detecting the existence and identifying the groups of data that are driven by the global common factors. The power of the test, together with the accuracy of various estimators, are demonstrated in extensive simulation settings. Estimation performance is also compared to other state-of-the-art methods for matrix time series by simulations and a set of U.S. equity data. All numerical works are done using the R package MMEFM, which is available at CRAN.
Random features provide a standard approach to scaling kernel methods to large datasets by approximating the kernel with a finite-dimensional feature map. Existing works of random feature methods primarily focus on well-specified models and light-tailed noise conditions, leaving a gap in understanding their performance under more general conditions, such as misspecified models and heavy-tailed noise. In this paper, we provide a comprehensive and improved analysis of kernel ridge regression with random features (KRR-RF) under general conditions. In the light-tailed noise setting, we improve upon the worst-case uniform bounds of the feature map used in the previous works and lead to the minimax optimal learning rates with $2r+γ>α_0$, where $(r,γ,α_0)$ denote the regularity of the target function, the effective dimension of the kernel, and the embedding index of the RKHS, respectively. In the heavy-tailed noise setting, we leverage a recently developed Hilbert-valued Fuk-Nagaev inequality and extend the same optimal results under finite weak moment conditions. The additional cost to achieve robustness is a mildly stronger explicit sample-size requirement than in the light-tailed case. For both settings, our results recover the conditions random features budgets (up to log factors) required in the previous works. In the most benign case when $α_0=γ$, KRR-RF are minimax optimal for all $r \in (0,1]$. To the best of our knowledge, this is the first work to establish optimal learning rates for KRR-RF under weak moment conditions in the misspecified setting. Numerical experiments and a real data analysis validate the derived learning rates and demonstrate their practical relevance in the presence of mis-specification and heavy-tailed noise.
Stochastic differential equation (SDE) models can reproduce trajectories while missing temporal dependence. We compare residual-based scores from a held-out path with scores from paths simulated by a fitted implementation, keeping fitted parameters, simulation rules and the scoring rule fixed. We give conditions for Monte Carlo calibration and bound conditional rejection error by the gap between conditional score CDFs. Residuals use zero-noise centres and covariance approximations that track noise propagation over observation intervals; report residuals also include reporting variation. Residual weights affect sensitivity and the influence of fitting error. For a quadratic energy score in a linear example with singular true covariance, whitening with nonsingular fitted covariance matrices leaves a nonvanishing score-CDF gap as nonzero coefficient estimates approach the truth. A fixed positive-definite weighting matrix makes the gap vanish. In a controlled SIR experiment using reference paths simulated at fitted parameters, adding susceptible-state endpoints to training or downweighting a low-variance residual reduces null rejection. A separate comparison using reference paths simulated at the true parameter shows reduced sensitivity to added susceptible-coordinate noise after downweighting. An exploratory weekly-report application incorporates reporting and reconstruction into reference simulations, revealing excess correlation between squared residuals in Czechia and persistently negative residuals in Finland.
In many applications, covariates are naturally organized into groups, yet conventional group-selection methods often retain all variables within selected groups, failing to achieve within-group sparsity. To address this limitation, we propose an adaptive two-stage sparse-group selection procedure for high-dimensional quantile regression. The proposed method integrates convolution-smoothed quantile loss with $\ell_{2,0}$ and $\ell_0$ penalties to identify important groups and relevant variables within them, achieving sparsity at both levels. We establish theoretical guarantees for sure screening and exact support recovery at both stages, as well as consistent model-size selection using information criteria, with explicit bounds on failure probabilities. The proposed procedure is computationally efficient by restricting optimization to selected groups and active variables. Simulation studies and a breast cancer gene-expression analysis demonstrate its competitive performance in variable selection, estimation accuracy, and computational efficiency, while identifying biologically relevant genes and delivering accurate quantile predictions.
In the literature on multiple testing methods that control the familywise error rate, a fundamental concept is "coherence". A procedure is coherent if it is "logically consistent" (Hochberg and Tamhane, 1987) in the sense that if a hypothesis is rejected, then all hypotheses that imply it are also rejected. If a method is not coherent, then it can be made coherent without increasing the familywise error rate. Recently, the concept of coherence was generalized to a much broader class of multiple testing procedures, namely all methods that provide simultaneous confidence bounds for false discovery proportions. Such procedures were defined to be coherent if they cannot be further improved through a process called "interpolation". However, we show that such a coherent procedure can sometimes be further improved by exploiting remaining logical inconsistencies. This leads to what we call a "perfectly coherent" method. Thus, interpolation does not necessarily lead to perfect coherence. However, we provide conditions under which it does. Fortunately, many methods based on interpolation turn out to be perfectly coherent.
Hierarchical time series forecasting requires generating accurate predictions that are coherent across multiple levels of aggregation. Existing methods either rely on two-stage pipelines with post-hoc reconciliation or on end-to-end models that often lack statistical guarantees. We introduce a unified framework that formulates hierarchical forecasting as a multivariate regression problem in a vector-valued reproducing kernel Hilbert space. By designing structured operator-valued kernels, our approach enforces coherence directly within the hypothesis space, ensuring that all predictions satisfy hierarchical constraints by construction. We further propose tailored kernel constructions for hierarchical data and derive excess risk bounds that characterize generalization performance. Empirical results on synthetic and real-world datasets demonstrate the relevance of the proposed method.
Randomly initialised neural networks of width $m$ converge to Gaussian processes as $m \to \infty$, with two types of quantitative bounds characterising this convergence. Finite-sample bounds compare the network with its Gaussian limit at a fixed finite sample set, on an arbitrary domain, in the quadratic Wasserstein distance $\mathcal W_2$, achieving a convergence rate of $m^{-1/2}$ with explicit constants. However, these constants grow with the number of samples, so the bound becomes uninformative as the sample set is refined. Functional bounds remove that dependence by comparing the two as random functions, at the price of spherical input, metrics weaker than $\mathcal W_2$, or slower rates. In this paper, we therefore lift a finite-sample bound to a functional one. The lifting refines the sample set until it is dense in the domain, so the bound must first be made independent of the number of samples. We achieve this by measuring the error per sample, which makes the constant uniform. For a polynomial activation and any compact input set, this gives a rate of $m^{-1/2}$ in the supremum norm, with an explicit constant.
This article proposes rational quadratic kriging, a kriging estimator that addresses the limitations of existing approaches: the mean reversion issue of ordinary kriging and the diverging variance issue of rational kriging. Through the rational quadratic kernel, we establish a connection between rational quadratic kriging and kernel regression. The connection allows us to prove desirable theoretical properties of the new estimator and further propose a novel estimation procedure that utilizes factor weights from kernel regression. We also derive expressions to enable efficient leave-one-out estimation and conformal prediction. Extensive simulations across both deterministic functions and noisy datasets show significant improvements in both prediction and uncertainty quantification. A fast implementation is available in the open-source R package RQK.
In this paper, we study the estimation of the signal subspace dimension $k$ in a $p$-variate signal-plus-noise complex Gaussian model, when the noise covariance matrix is assumed to be diagonal or Toeplitz. We consider the high-dimensional asymptotic regime in which both the data dimension $p$ and the available sample size $n$ diverge to infinity at the same rate, while we also allow $k$ to diverge at a slower speed. In this context, the classical estimates of the literature, consistent in the high-dimensional regime, often rely on the a priori knowledge of an upper bound $k_*$ of $k$, which can be misspecified in practice. Our proposed approach relies on building estimates of the tracial moments of the noise covariance matrix using the $p-k_*$ smallest eigenvalues of the sample covariance matrix (SCM), which are in turn used to build a threshold above which the $k$ largest eigenvalues are isolated. The resulting estimate of $k$ is proved to be consistent in the high-dimensional regime, even for the misspecified situations where $k_* < k$. Numerical experiments in different challenging scenarios illustrate the improvement of the proposed approach over the aforementioned methods.
We study stochastic bandits without knowledge of reward scales or tail parameters. Almost-sure regret guarantees describe eventual learning but do not quantify expected performance or the risk of large regret in a finite run. We introduce $\texttt{EL-F}$, an anytime empirical-likelihood (EL) extension of the $g$-Forcing policy of Cowan and Katehakis (2020), and characterize its performance at three levels. First, for integrable rewards and non-decreasing, diverging exploration schedules $g$ with vanishing increments, $\texttt{EL-F}$ preserves the $g$-Forcing asymptotics for pseudo-regret $R(T)$: $R(T)/g(T)\to\sum_{i:Δ_i>0}Δ_i$ almost surely as $T\to\infty$, where $Δ_i$ is the suboptimality gap of arm $i$. Second, on every sub-exponential instance, for non-decreasing $g(t)=ω(\log t)$ satisfying a mild feasibility condition, expected pseudo-regret is at most $g(T)\sum_{i:Δ_i>0}Δ_i+C_{\underline{P}}$, where $C_{\underline{P}}$ is an instance-dependent constant independent of $T$. Moreover, we also establish finite-time tail bounds with super-polynomial decay in regret beyond the forced-exploration budget. For slowly varying $g$, we determine the exact exponent of linear regret deviations on two-armed Gaussian instances with common variance and establish a matching lower bound among policies satisfying the same expected-regret budget on every such Gaussian instance. These quantitative guarantees also cover the greedy variant obtained at zero EL radius. A fixed positive radius introduces local optimism while preserving the leading asymptotic regret. We establish coupled pathwise dominance over the greedy variant in two-armed bandits with a deterministic suboptimal arm. Experiments demonstrate reductions in prolonged suboptimal allocation compared with $g$-Forcing on challenging instances.
Nested dichotomies (NDs) reduce multiclass classification to a tree of binary classifiers. Tree structure can change accuracy, bias, and variance, but the full space of structures has not been mapped in aggregate. We introduce ndscape, a two-dimensional mapping of ND structures based on depth-weighted structural tokens. We use this embedding to conduct what is, to our knowledge, the first bias-variance analysis across the ND structure space. In simulations with 6 to 10 classes, NDs with similar tree structures exhibit similar squared bias, variance, and squared probability loss, regardless of whether the binary classifier is logistic regression or a decision tree. Variance, however, follows mean root-to-leaf depth when logistic regression is used at each node, so flatter, more balanced NDs are more stable, whereas with decision-tree node models, variance is not associated with depth and differs little across ND structures. On six datasets with 6 to 19 classes, held-out accuracy, log loss, and prediction variance remain spatially organised on the same topology-only map. None of the six standard construction heuristics tested selects a tree from the best-performing cluster on any dataset, although some select high-performing individual trees. The ndscape embedding turns ND selection from one opaque construction step into a readable comparison of class groupings, exposing the upper-level splits that distinguish high-performing and stable regions.