337 papers
Disentangled Deep Priors for Bayesian Inverse Problems
We propose a structured prior for high-dimensional Bayesian inverse problems based on a disentangled deep generative model whose latent space is partitioned into auxiliary variables aligned with known and interpretable physical parameters and residual variables capturing remaining unknown variability. This yields a hierarchical prior in which interpretable coordinates carry domain-relevant uncertainty while the residual coordinates retain the flexibility of deep generative models. By linearizing the generator, we characterize the induced prior covariance and derive conditions under which the posterior exhibits approximate block-diagonal structure in the latent variables, clarifying when representation-level disentanglement translates into a separation of uncertainty in the inverse problem. We formulate the resulting latent-space inverse problem and solve it using MAP estimation and Markov chain Monte Carlo (MCMC) sampling. On elliptic PDE inverse problems, such as conductivity identification and source identification, the approach matches an oracle Gaussian process prior under correct specification and provides substantial improvement under prior misspecification, while recovering interpretable physical parameters and producing spatially calibrated uncertainty estimates.
2604.02304Apr 2026
View
Simulation and Analysis of Multifractional Stochastic Processes with R Package Rmfrac
Brownian motion and fractional Brownian motion have been widely applied in statistical modeling in finance, telecommunication, network traffic, neuroscience, physics, and other fields. More realistic models for real time series data, such as multifractional processes, generalize these classical models by allowing their regularity to vary over time. A new class of Gaussian Haar-based multifractional processes, which utilizes the Haar wavelet series representation, was recently introduced. It significantly extends the range of available models by incorporating more general classes of Hurst functions. The Rmfrac package was developed to simulate multifractional time series. The package also comprises several functions for the analysis and visualization of time series. It includes the estimation of the Hurst function and local fractal dimension, clustering realizations and computing various geometric statistics of these time series. The package also offers a Shiny application to visualize simulation and estimation results. The article presents an overview of the Rmfrac package and exemplifies its main functionalities.
2604.01580Apr 2026
View
Quantum Statistical Bootstrap
The bootstrap is a foundational tool in statistical inference, but its classical implementation relies on Monte Carlo resampling, introducing approximation error and incurring high computational cost -- especially for large datasets and complex models. We present the Quantum Bootstrap (QBOOT), a quantum algorithm that computes the ideal bootstrap estimate exactly by encoding all possible resamples in quantum superposition, evaluating the target statistic in parallel, and extracting the aggregate via quantum amplitude estimation. Under mild circuit efficiency assumptions, QBOOT achieves a near-quadratic speedup over the classical bootstrap in approximating the ideal estimator, independent of the statistic or underlying distribution. We provide a rigorous theoretical analysis of its statistical error properties -- addressing a gap in the quantum algorithms literature -- and validate our results through experiments on the IBM quantum simulator for the sample mean problem. Our findings demonstrate that QBOOT preserves the asymptotic properties of the ideal bootstrap while substantially improving computational efficiency and precision, establishing a scalable and principled framework for quantum statistical inference.
2604.00951Apr 2026
View
Implementation and Workflows for INLA-Based Approximate Bayesian Structural Equation Modelling
Bayesian structural equation modelling (BSEM) offers many advantages such as principled uncertainty quantification, small-sample regularisation, and flexible model specification. However, the Markov chain Monte Carlo (MCMC) methods on which it relies are computationally prohibitive for the iterative cycle of specification, criticism, and refinement that careful psychometric practice demands. We present INLAvaan, an R package for fast, approximate Bayesian SEM built around the Integrated Nested Laplace Approximation (INLA) framework for structural equation models developed by Jamil & Rue (2026, arXiv:2603.25690 [stat.ME]). This paper serves as a companion manuscript that describes the architectural decisions and computational strategies underlying the package. Two substantive applications -- a 256-parameter bifactor circumplex model and a multilevel mediation model with full-information missing-data handling -- demonstrate the approach on specifications where MCMC would require hours of run time and careful convergence work. In constrast, INLAvaan delivers calibrated posterior summaries in seconds.
2604.00671Apr 2026
View
2603.29724A new gradient-free active subspace estimation method with application to rare event probability estimation
To reduce the cost of estimating the probability of a rare event involving a very large number of random parameters, we propose a new strategy for dimension reduction coupled with a surrogate model for the expensive part of the algorithm. To this end, we extend the Ordinary Kriging Active Subspace (OK-AS) method into a sequential version. Our approach consists of iteratively re-estimating the active subspace using a Kriging surrogate trained in a rotated coordinate system until the active subspace stabilises. This method allows for a reduction in prediction error and a better approximation of the active subspace on a benchmark of test problems. Furthermore, we integrate our algorithm into an efficient pre-existing approach for estimating the probability of a rare event. This approach is based on learning the active subspace associated with the random event whose probability is to be estimated. The sequential learning of an importance sampling density is necessary and corresponds to the expensive part of this strategy. To circumvent this issue, we integrate our sequential OK-AS version into the estimation of the importance sampling density. The numerical results indicate that our method allows for reducing the cost required to obtain a precise estimate of the rare event probability.
2603.29724Mar 2026
View
A Hybrid NUTS-Gibbs Sampler with State Space Marginalization for Estimation of Dynamic Structural Equation Models with Binomial Outcomes
Dynamic structural equation modeling (DSEM) is widely used for analyzing intensive longitudinal data (ILD). Although many ILD have categorical (Bernoulli or binomially distributed) responses, currently available Metropolis-within-Gibbs samplers for estimating DSEMs are limited to using the probit link and the Bernoulli distribution. These samplers scale poorly with increasing model complexity and/or data size. Here, we present a hybrid sampler -- alternating between one step of the No-U-Turn Sampler (NUTS) and one Gibbs step -- which solves both of these problems: the Gibbs step naturally handles Pólya-Gamma distributed latent variables arising from binomially distributed responses with a logit link, and the NUTS step utilizes a Kalman filter to exactly marginalize over latent states, alleviating the need to sample these variables. We demonstrate in simulation experiments that the proposed sampler is more efficient than alternative algorithms, and that it makes DSEM estimation with binomial data feasible for larger data and models than what has previously been possible. We also illustrate its use in an example application of predicting panic attacks.
2603.29647Mar 2026
View
Retrospective Economic Evaluation of Group Testing in the COVID-19 Pandemic
Surveillance of diseases in a pandemic is an important part of public health policy. Diagnostic testing at the individual level is often infeasible due to resource constraints. To circumvent these constraints, group testing can be applied. The economic cost evaluation from the payer's perspective typically focuses only on deterministic costs which overlooks the substantial economic impact of productivity losses resulting from quarantine and workplace disruptions. The objective of this article is to develop a mathematical model for a retrospective economic evaluation of group testing that incorporates both deterministic costs and income-based economic loss. Group testing algorithms are revisited and simulated at optimized pool sizes to determine the required number of tests. Income data from the German Socio-Economic Panel are integrated into a mathematical model to capture the economic loss. Afterward, hybrid Monte Carlo experiments are conducted by evaluating the economic cost in the Coronavirus disease 2019 pandemic in Germany. Monte Carlo experiments show that the optimal choice of group testing algorithms changes substantially when income-based economic losses are included. Evaluations considering only deterministic costs systematically underestimate the total economic cost. Algorithms with a longer quarantine duration are less attractive than shorter quarantine duration if income-based economic loss is accounted for. The findings show that current evaluations underestimate the true economic cost. Group testing algorithms with shorter duration and fewer stages are preferred, even when they require a larger number of tests. These results underscore the importance of incorporating income-based economic loss into a mathematical model.
2603.28930Mar 2026
View
Generalized and Scalable Deep Gaussian Process Emulation
Gaussian process (GP) emulators have become essential tools for approximating complex simulators, significantly reducing computational demands in optimization, sensitivity analysis, and model calibration. While traditional GP emulators effectively model continuous and Gaussian-distributed simulator outputs with homogeneous variability, they typically struggle with discrete, heteroskedastic Gaussian, or non-Gaussian data, limiting their applicability to increasingly common stochastic simulators. In this work, we introduce a scalable Generalized Deep Gaussian Process (GDGP) emulation framework designed to accommodate simulators with heteroskedastic Gaussian outputs and a wide range of non-Gaussian response distributions, including Poisson, negative binomial, and categorical distributions. The GDGP framework leverages the expressiveness of DGPs and extends them to latent GP structures, enabling it to capture the complex, non-stationary behavior inherent in many simulators while also modeling non-Gaussian simulator outputs. We make GDGP scalable by incorporating the Vecchia approximation for settings with a large number of input locations, while also developing efficient inference procedures for handling large numbers of replicates. In particular, we present methodological developments that further enhance the computation of the approach for heteroskedastic Gaussian responses. We demonstrate through a series of synthetic and empirical examples that these extensions deliver the practical application of GDGP emulators and a unified methodology capable of addressing diverse modeling challenges. The proposed GDGP framework is implemented in the open-source R package dgpsi.
2603.24538Mar 2026
View
Uncertainty Quantification of Spline Predictors on Compact Riemannian Manifolds
To predict smooth physical phenomena from observations, spline interpolation provides an interpretable framework by minimizing an energy functional associated with the Laplacian operator. This work proposes a methodology to construct a spline predictor on a compact Riemannian manifold, while quantifying the uncertainty inherent in the classical deterministic solution. Our approach leverages the equivalence between spline interpolation and universal kriging with a specific covariance kernel. By adopting a Gaussian random field framework, we generate stochastic simulations that reflect prediction uncertainty. However, on compact manifolds, the covariance kernel depends on the generally unknown spectrum of the Laplace-Beltrami operator. To address this, we introduce a finite element approximation based on a triangulation of the manifold. This leads to the use of intrinsic Gaussian Markov Random Fields (GMRF) and allows for the incorporation of anisotropies through local modifications of the Riemannian metric. The method is validated using a temperature study on a sphere, where the operator's spectrum is known, and is further extended to a test case on a cylindrical surface.
2603.24259Mar 2026
View
Adaptive Gaussian Process Search for Simulation-Based Sample Size Estimation in Clinical Prediction Models: Validation of the pmsims R Package
Background: Determining an adequate sample size is essential for developing reliable and generalisable clinical prediction models, yet practical guidance on selecting appropriate methods remains limited. Existing analytical and simulation-based approaches often rely on restrictive assumptions and focus on mean-based criteria. We present and validate pmsims, an R package that uses Gaussian process surrogate modelling to provide a flexible and computationally efficient simulation-based framework for sample size determination across diverse prediction settings. Methods: We conducted a comprehensive simulation study with two aims. First, we compared three search engines implemented in pmsims: a Gaussian process-based adaptive method, a deterministic bisection method, and a hybrid approach, across binary, continuous, and survival outcomes. Second, we benchmarked the best-performing pmsims engine against existing analytical (pmsampsize) and simulation-based (samplesizedev) methods, evaluating recommended sample sizes, computational time, and achieved performance on large independent validation datasets. Results: The Gaussian process-based method consistently produced the most stable sample size estimates, particularly in low-signal, high-dimensional settings. In benchmarking, pmsims achieved performance close to prespecified targets across all outcome types, matching simulation-based approaches and outperforming analytical methods in more challenging scenarios. Conclusions: pmsims provides an efficient and flexible framework for principled sample size planning in clinical prediction modelling, requiring fewer model evaluations than non-adaptive simulation approaches.
2603.23688Mar 2026
View
Neural Inference Functions for Margins for Time Series Copula Models
Copula models are widely employed in multivariate time series analysis because they permit flexible modelling of marginal distributions independently of the dependence structure, which is fully characterised by the copula function. However, Bayesian inference with these models becomes computationally demanding as the number of variables in the time series increases. Motivated by the classical inference functions for margins (IFM) approach, we propose a new neural-network based inference framework for estimating parameters in copula models, termed the neural inference functions for margins (N-IFM). N-IFM enables rapid parameter estimation for new data, fast sequential prediction, and efficient model comparison via time-series validation. We assess the performance of N-IFM using both simulated and real datasets and compare it to Hamiltonian Monte Carlo, demonstrating substantial computational gains with comparable inferential accuracy.
2603.21075Mar 2026
View
A practical introduction to ODE modelling in Stan for biological systems
Integrating dynamical systems models with time series data is a central part of contemporary mathematical biology. With the rich variety of available models and data, numerous methods and computational tools have been developed for these purposes. One such tool is Stan, a freely available and open-source probabilistic programming framework that provides efficient methods for estimating model parameters from data using computational Bayesian inference algorithms. Stan includes built-in mechanisms for working with ordinary differential equation (ODE) models, which are widely used in mathematical biology and related fields to study simulated, experimental, and real-world systems that change over time. Through step-by-step worked examples, including both pedagogical toy models and applications with real data, this article provides a practical, self-contained introduction to performing parameter estimation and model evaluation for first-order linear and nonlinear ODE models in Stan. The article also explains key statistical methods that underpin Stan and discusses computational Bayesian modelling in the context of biological applications.
2603.20343Mar 2026
View
Adaptive Nonlinear Data Assimilation through P-Spline Triangular Measure Transport
Non-Gaussian statistics are a challenge for data assimilation. Linear methods oversimplify the problem, yet fully nonlinear methods are often too expensive to use in practice. The best solution usually lies between these extremes. Triangular measure transport offers a flexible framework for nonlinear data assimilation. Its success, however, depends on how the map is parametrized. Too much flexibility leads to overfitting; too little misses important structure. To address this balance, we develop an adaptation algorithm that selects a parsimonious parametrization automatically. Our method uses P-spline basis functions and an information criterion as a continuous measure of model complexity. This formulation enables gradient descent and allows efficient, fine-scale adaptation in high-dimensional settings. The resulting algorithm requires no hyperparameter tuning. It adjusts the transport map to the appropriate level of complexity based on the system statistics and ensemble size. We demonstrate its performance in nonlinear, non-Gaussian problems, including a high-dimensional distributed groundwater model.
2603.19058Mar 2026
View
Preconditioning Hamiltonian Monte Carlo by minimizing Fisher Divergence
Although Hamiltonian Monte Carlo (HMC) scales as O(d^(1/4)) in dimension, there is a large constant factor determined by the curvature of the target density. This constant factor can be reduced in most cases through preconditioning, the state of the art for which uses diagonal or dense penalized maximum likelihood estimation of (co)variance based on a sample of warmup draws. These estimates converge slowly in the diagonal case and scale poorly when expanded to the dense case. We propose a more effective estimator based on minimizing the sample Fisher divergence from a linearly transformed density to a standard normal distribution. We present this estimator in three forms, (a) diagonal, (b) dense, and (c) low-rank plus diagonal. Using a collection of 114 models from posteriordb, we demonstrate that the diagonal minimizer of Fisher divergence outperforms the industry-standard variance-based diagonal estimators used by Stan and PyMC by a median factor of 1.3. The low-rank plus diagonal minimizer of the Fisher divergence outperforms Stan and PyMC's diagonal estimators by a median factor of 4.
2603.18845Mar 2026
View
Bridging Theory and Practice in Efficient Gaussian Process-Based Statistical Modeling for Large Datasets
Geostatistics is a branch of statistics concerned with stochastic processes over continuous domains, with Gaussian processes (GPs) providing a flexible and principled modelling framework. However, the high computational cost of simulating or computing likelihoods with GPs limits their scalability to large datasets. This paper introduces the piecewise continuous Gaussian process (PCGP), a new process that retains the rich probabilistic structure of traditional GPs while offering substantial computational efficiency. As will be shown and discussed, existing scalable approaches that define stochastic processes on continuous domains -- such as the nearest-neighbour GP (NNGP) and the radial-neighbour GP (RNGP) -- rely on conditional independence structures that effectively constrain the measurable space on which the processes are defined, which may induce undesirable probabilistic behaviour and compromise their practical applicability, particularly in complex latent GP models. The PCGP mitigates these limitations and provides a theoretically grounded and computationally efficient alternative, as demonstrated through numerical illustrations.
2603.18324Mar 2026
View
Sequential Bayesian Experimental Design for Prediction in Physical Experiments Informed by Computer Models
In many scientific and engineering domains, physical experiments are often costly, non-replicable, or time-consuming. The Kennedy and O'Hagan (KOH) model framework has become a widely used approach for combining simulator runs with limited experimental observations. Under a Bayesian implementation, the simulator output, model discrepancy, and observation noise are jointly modeled by coupled Gaussian processes, followed by coherent posterior inference and uncertainty quantification. This work presents a genuinely sequential Bayesian experimental design (BED) framework explicitly aimed at improving the predictive performance of the KOH model. We employ a mutual information (MI)-based criterion and develop a hybrid variant that integrates it with measures of local model complexity, leading to significantly more efficient design decisions. We further show theoretically that the MI-based criterion is more comprehensive and robust than the classical integrated mean squared prediction error (IMSPE) minimization criterion, especially when the model is highly uncertain in the early stages of the experiment. To mitigate the computational burden of fully Bayesian inference and the ensuing BED process, we propose two acceleration strategies - Gaussian Mixture Compression and Schur complement and rank-one update - which together substantially reduce runtime. Finally, we demonstrate the effectiveness of the proposed methods through both a synthetic example and a real biochemical case study, and compare them against several classical design criteria under sequential (offline) and adaptive (online) BED settings.
2603.16756Mar 2026
View
Fast Multitask Gaussian Process Regression
Gaussian process (GP) regression is a powerful probabilistic modeling technique with built-in uncertainty quantification. When one has access to multiple correlated simulations (tasks), it is common to fit a multitask GP (MTGP) surrogate which is capable of capturing both inter-task and intra-task correlations. However, with a total of $N$ evaluations across all tasks, fitting an MTGP is often infeasible due to the $\mathcal{O}(N^2)$ storage and $\mathcal{O}(N^3)$ computations required to store, solve a linear system in, and compute the determinant of the $N \times N$ Gram matrix of pairwise kernel evaluations. In the single-task setting, one may reduce the required storage to $\mathcal{O}(N)$ and computations to $\mathcal{O}(N \log N)$ by fitting "fast GPs" which pair low-discrepancy design points from quasi-Monte Carlo to special kernel forms which yields nicely structured Gram matrices, e.g., circulant matrices. This article generalizes fast GPs to fast MTGPs which pair low-discrepancy design points for each task to special product kernel forms which yields nicely structured block Gram matrices, e.g., circulant block matrices. An algorithm is presented to efficiently store, invert, and compute the determinant of such Gram matrices with optionally different sampling nodes and different sample sizes for each task. Derivations for fast MTGP Bayesian cubature are also provided. A GPU-compatible, open-source Python implementation is made available in the FastGPs package (https://alegresor.github.io/fastgps/). We validate the efficiency of our algorithm and implementation compared to standard techniques on a range of problems with low numbers of tasks and large sample sizes.
2603.16014Mar 2026
View
SEMMS with Random Effects: A Mixed-Model Extension for Variable Selection in Clustered and Longitudinal Data
SEMMS (Scalable Empirical-Bayes Model for Marker Selection) is a variable-selection procedure for generalized linear models that uses a three-component normal mixture prior on regression coefficients. In its original form, SEMMS assumes that all observations are independent. Many real-world datasets, however, arise from repeated-measures or clustered designs in which observations within the same subject are correlated. Ignoring this correlation inflates the apparent residual variance and can severely degrade variable-selection performance. We extend SEMMS to accommodate random intercepts, random slopes, or both, via an alternating coordinate-ascent algorithm. After each round of fixed-effect variable selection, the subject-level best linear unbiased predictors (BLUPs) are updated with \texttt{lmer} (Gaussian) or \texttt{glmer} (non-Gaussian); the fixed-effect step then operates on the random-effect-adjusted response. We describe the algorithm, evaluate its performance in three Gaussian simulation studies spanning a range of signal strengths, random-effect magnitudes, and sample/predictor-space regimes, and present a semi-synthetic real-data example. We further extend the framework to non-Gaussian families (Poisson, binomial) via an IRLS working-response adaptation: at each outer iteration the fixed-effects step uses the RE-adjusted working response computed from the current \texttt{glmer} fitted values rather than the raw response. When the fixed-effect signal is strong relative to the random-effect variance, both the original and extended procedures perform comparably. When the random-effect variance dominates -- the scenario most likely to cause plain SEMMS to fail -- the mixed-model extension recovers the exact true predictor set in 93\% of simulated datasets (Gaussian), 61\% (Poisson), and 65\% (binomial), compared with 1\%, 45\%, and 39\% for plain SEMMS respectively.
2603.15902Mar 2026
View2603.13622The Continuous Rank Probability Score of a Generalized Beta-Prime Distribution and Some Special Cases
This working paper describes new results in derivations of the Continuous Ranked Probability Score of a generalized beta-prime distribution and several special cases, such as the Dagum distribution and Singh-Maddala distribution. Comparison with Monte Carlo estimates is also presented.
2603.13622Mar 2026
View
Edgewise Envelopes Between Balanced Forman and Ollivier-Ricci Curvature
Evaluating Ollivier-Ricci (OR) curvature on large-scale graphs is computationally prohibitive due to the necessity of solving an optimal transport problem for every edge. We bypass this computational bottleneck by deriving explicit, two-sided, piecewise-affine transfer moduli between the transport-based OR curvature and the combinatorial Balanced Forman (BF) curvature introduced by Topping et al. By constructing a lazy transport envelope and augmenting the Jost and Liu bound with a cross-edge matching statistic, we establish deterministic bounds for $\mathfrak{c}_{OR}(i,j)$ parameterized by 2-hop local graph combinatorics. This formulation reduces the edgewise evaluation complexity from an optimal transport linear program to a worst-case $\mathcal{O}(\max_{v \in V} \operatorname{deg}(v)^{1.5})$ time, entirely eliminating the reliance on global solvers. We validate these bounds via distributional analyses on canonical random graphs and empirical networks; the derived analytical bands enclose the empirical distributions independent of degree heterogeneity, geometry, or clustering, providing a scalable, computationally efficient framework for statistical network analysis.
2603.13535Mar 2026
ViewPage 1 of 17
