arXane › Statistics › Statistics
Wednesday, 4 November 2026: 2 new papers, 7 cross-lists and 6 replacements in Computation (stat.CO), in arXiv's order, with full abstracts. The same papers as arXiv's own listing, set as a newspaper; every title links to arXiv.
Computing phylogenetic likelihoods becomes challenging when the continuous-time Markov chain has a large state space. At each branch, pruning algorithms require terms of the form $e^{Q t}\cdot B$, where $Q$ is an $S \times S$ rate matrix (e.g., for $S$ states), $t$ is a branch length, and the columns of matrix $B$ (e.g., $B \in \mathbb{R}_{+}^{S \times L}$) are partial-likelihood vectors for $L$ site patterns. Many implementations calculate the dense transition matrix $P(t)=e^{Q t}$ via eigendecomposition and multiply it by $B$. This strategy is efficient for small state spaces, but it can be inefficient for large sparse rate matrices because forming $P(t)$ or its eigen vector matrices destroys sparsity and its cost may not be well amortized when $L$ is small. We implement an alternative likelihood kernel based on the Al-Mohy and Higham (2011) algorithm for computing the action of the matrix exponential directly, avoiding explicit construction of $P(t)$ or its eigen vector matrices. The resulting method takes full advantage of the sparsity of $Q$. We implement this method as a CPU-based extension to the BEAGLE library and evaluate it on two large-state molecular models: a Markov-modulated codon model of heterotachy with $M$ selective regimes and $61\cdot M$ possible states, and a codon-pair model of interlocus gene conversion between two paralogs with $61^2=3\,721$ possible states. Our results characterize when action-based likelihood computation improves performance over calculating $e^{Q t}$ directly. We show that this method enables fast, scalable inference for complex models, removing a key computational barrier to the study of evolution with large state space models.
Inference in stochastic reaction-network models, such as the susceptible-exposed-infectious-recovered (SEIR) epidemic model, is crucial for understanding the dynamics of interacting systems in epidemiology, ecology, and systems biology. These models are typically represented as Markov jump processes (MJPs), with discrete, noisy observations and intractable likelihoods. Particle Markov chain Monte Carlo (PMCMC) offers a natural mechanism for inference on these hidden Markov models, but, even in moderate dimensional state spaces, particle degeneracy typically leads to poor mixing of the PMCMC algorithm. This paper focuses on improving the efficiency of particle Gibbs for inference on reaction networks by addressing the degeneracy problem. Building on recent work on the exchangeable particle Gibbs (xPG) sampler for diffusions, this paper develops a novel version of xPG tailored to discrete-state-space reaction networks, where randomness is driven by Poisson processes rather than Brownian motion. We focus on the $\tau$-leap discretisation of the MJPs. The proposed method retains the exchangeability framework of xPG while adapting it to the structural and computational challenges specific to reaction networks. The new method is applied to several examples, including four coupled SEIR models, where it substantially improves mixing efficiency compared with standard particle Gibbs. In this example, xPG improves sample efficiency by more than an order of magnitude relative to standard particle Gibbs, making inference feasible in hours rather than weeks. As part of our tuning strategy, we provide a mechanism for tuning the number of particles that is applicable more generally, to any particle Gibbs algorithm.
We study the multi-modal Laplace approximation, given by a weighted Gaussian mixture centred at the local minima of the negative log-density, for sampling from probability measures that concentrate around finitely many modes as an inverse-temperature parameter $n$ grows. Such distributions occur in Bayesian inference in the small noise or large data regime. We derive explicit, non-asymptotic bounds on the Hellinger distance between the multi-modal Laplace approximation and the concentrating posterior. The bound decays like $n^{-1/2}$ with $n$ and the constants depend only polynomially on the dimension $d$. Moreover, we analyse the Laplace mixture as a proposal distribution for self-normalised importance sampling and independent Metropolis--Hastings (IMH) targeting concentration multimodal distributions. We derive asymptotic results, as $n\to\infty$, on the effective sample size in case of importance sampling as well as on the average acceptance probability, lag-one autocorrelation, and spectral gap in case of IMH. Numerical experiments illustrate and verify our theoretical findings.
Distributed statistical learning has recently attracted significant attention owing to its computational efficiency and scalability in large-scale systems. However, the non-smooth quantile loss poses substantial challenges to the development of efficient decentralized algorithms for quantile regression. Existing methods may either suffer from slow convergence or rely on smoothing-based approximations that require additional stringent assumptions, particularly in high-dimensional settings. To address these issues, we propose a decentralized framework that integrates gradient tracking with iterative hard thresholding for sparse quantile estimation. The proposed method avoids smoothing approximations and relaxes the initialization requirement. Theoretically, we show that the proposed algorithm achieves linear convergence under stated regularity conditions and that the resulting estimator attains the same order of statistical accuracy as its centralized counterpart. We further develop a decentralized information criterion for selecting the working sparsity level and establish support screening guarantees. Numerical experiments and a real-data study illustrate the finite-sample performance of the proposed method.
After variable selection, uncertainty quantification often treats the selected model as fixed, although it may change under repeated sampling. We investigate whether this variability can be incorporated by repeating Cox-Lasso selection within a model-based bootstrap after unpenalised Cox refitting. We compare conventional post-selection Cox Wald intervals and a model-based bootstrap without re-selection, both of which keep the selected predictor set fixed, with two re-selection procedures: one retains the observed tuning-parameter value, whereas the other reapplies the cross-validation minimum (CV-min) tuning rule within each bootstrap sample. Against the data-generating coefficients, procedures holding the selected predictor set fixed frequently showed substantial undercoverage. Re-selection generally moved coverage closer to the nominal level, although performance varied across simulation regimes and coefficients. Sensitivity analyses evaluating coverage against the population Cox coefficients corresponding to each originally selected predictor set gave similar aggregate coverage for the primary interval comparison. In the NKI70 breast-cancer data, the four selected genes appeared together in only 0.6% of bootstrap samples under re-selection at the observed tuning-parameter value. These results suggest that selection instability can be an important component of post-selection coefficient uncertainty and that bootstrap re-selection can help incorporate this variability.
Risk and reliability analysis are fundamental to the design, assessment, and operation of safety-critical complex systems; however, their integration into next-generation computational tools is constrained by the cost of uncertainty quantification with expensive computational models. The computational challenge is particularly prohibitive for tail-sensitive quantities, such as probabilities of failure, conditional values-at-risk, and buffered probabilities of failure, where accurate estimation can require a large number of model evaluations. Multifidelity methods provide a framework for addressing this computational bottleneck by combining information from models of varying accuracy and cost while retaining the predictive accuracy of the high-fidelity model. This review summarizes the state of the art in methods that enable scalable risk and reliability analysis using multifidelity control variates, multifidelity importance sampling, and multifidelity active learning.
The Ising model, a type of Markov Random Field, from statistical physics has a long history of modeling spatial and network data as well as particle interactions. In this work, we consider inserting the Ising model into the framework of a Hidden Markov Model, a type of Bayesian Network, to model temporal state changes in a sequence of Ising observations. While this is entirely a discrete probability model, a main challenge is the fitting of such a probabilistic graphical model due to the exponential intractability of normalizing constants within the Ising model. To circumvent this problem, we consider the usage of maximum pseudolikelihood estimators that have a long history in Markov Random Fields. Validity of the pseudolikelihood Baum-Welch algorithm, model identifiability under likelihood, and consistency for large networks is also proved. The proposed Hidden Markov Ising Model is further extended to incorporate exogenous variables resulting in a hybrid hidden Markov autologistic model with graphical dependencies integrated into the model. This model is applied to a five year sequence of weekly gun violence events across the 3109 counties of the continental United States to identify state changes in gun violence activity from 2021 to 2025.
Importance sampling is a building block of modern statistical computing, but the weighted estimate is prone to high or infinite variance. Moreover, its empirical standard error is often unreliable, leaving importance sampling without practical uncertainty quantification. Motivated by multilevel regression and poststratification (MRP) in survey sampling, we propose post-stratified importance sampling (PIS) for estimating target expectations from proposal draws. PIS reformulates arbitrarily high-dimensional weighting into two one-dimensional probabilistic modeling problems: regression for the conditional mean given the importance ratio and marginal density estimation for the importance ratio. Different modeling choices yield a family of PIS estimators, encompassing ordinary importance sampling as a special case. We develop a model-based implementation that combines one-dimensional Gaussian process regression for the conditional mean and score-informed heavy-tailed density estimation for the marginal law, using both the observed weights and the available density evaluations. We study the oracle gains of PIS and establish consistency under general and model-specific conditions. Simulations and applications to cross-validation demonstrate improvements in estimation accuracy and uncertainty quantification.
Implicit generative models based on deep neural networks have revolutionized many scientific fields, from weather research to language generation. They successfully model complex random objects arising in real-world problems. In this work, we propose a sample-efficient training mechanism for these models. Our approach leverages regularized inference with shrinkage priors and generalized Bayes. In the large $p$ (parameter dimension) and small $n$ (sample size) setting, sparsity-inducing priors have been successfully used for linear regression when $p \gg n$. However, their application to implicit generative models remains challenging due to the intractability of the likelihood function. To address this, we adapt a scoring rule posterior based on a strictly proper scoring rule, following the generalized Bayes framework. We then propose a block SGMCMC within Gibbs sampling mechanism. This enables efficient exploration of high-dimensional parameter spaces by identifying a sparse representation of neural implicit generative models. We demonstrate the performance of our proposed method through experiments on $p \gg n$ linear regression and an application of neural generative models.