Papers

In sampling problems, gradient-based schemes such as Langevin Monte Carlo (LMC) mix faster than non-gradient-based methods, but their applicability is limited by access to the gradient of the target log-density. In practice, gradients are often unavailable and function evaluations are noisy, e.g., stochastic simulators or black-box simulators, so we propose LMC-SPSA with noise, which approximates the gradient of the target log-density using two noisy function evaluations per iteration. We prove, under noisy gradient estimates, that LMC-SPSA converges in distribution by proving the convergence in Wasserstein distance. Furthermore, we construct a diminishing step-size schedule that still drives the Wasserstein error bound to convergence, extending convergence guarantees beyond the constant-step setting. Further, we sharpen the dominant dimension dependence of the Wasserstein error from $O(p^4)$ to $O(p^2)$ (with $p$ denoting the dimension), and support this analysis with numerical results. We show that LMC-SPSA achieves $W_2$-accuracy $\varepsilon$ with total noisy-oracle complexity of $O(p/\varepsilon^2+\delta^2p^3/\varepsilon^3)$, where $\delta$ is the paired-noise level. This improves the noise-dependent accuracy scaling relative to the ZO-LMC method of Roy et al. We further establish asymptotically vanishing Wasserstein error as the number of iterations $\to\infty$ under diminishing step-size and perturbation sequences and derive an explicit convergence rate for a balanced schedule under noisy zeroth-order feedback. Empirical experiments are conducted to verify the performance of LMC-SPSA with noise. We provide an oracle-budget-matched comparison with the ZO-LMC method, showing smaller empirical sampling errors under the same function-evaluation budget.

This study investigated the behavior of Principal Component Analysis (PCA) when applied to datasets with extremely large numbers of observations. Although statistical theory suggests that sampling error diminishes and sample estimates converge toward their population values as sample size increases, relatively little empirical evidence exists regarding the behavior of PCA at scales measured in billions or trillions of observations. Three datasets were analyzed: a 10-billion observation random dataset, a 1-trillion observation random dataset, and a 10-billion observation engineered dataset designed to contain three latent factors. Results showed that the PCA solutions obtained from the 10BillionRandom and 1TrillionRandom datasets were nearly identical, indicating substantial stability of PCA at extremely large sample sizes. In contrast, the engineered dataset produced three dominant principal components that accounted for 99.996% of the total standardized variance and successfully recovered the intended latent-factor structure. These findings suggest that PCA solutions converge rapidly at very large sample sizes and suggest that PCA solutions may reach practical convergence well before sample sizes reach the trillions. These findings have implications for large-scale applications in fields such as remote sensing, digital mapping, environmental modeling, and other domains where datasets routinely contain millions or billions of observations.

We propose a new unit of analysis for longitudinal data: the Latent Memory Table. The scientific contribution is not the encoder. It is that table, treated as a reusable statistical object on the same footing as a matrix of principal-component scores, a table of estimated random effects, or a table of predicted probabilities. We estimate a statistical table that summarizes recent longitudinal history and is intended to be stored, queried, analysed and reused throughout the statistical workflow. A memory operator maps each masked windowed history to a finite-dimensional state; collecting those states with uncertainty yields the Latent Memory Table. Validation is organized around six properties—recoverability, personalization, temporal coherence, interpretability, stability and reusability—summarized by a composite quality index \(Q\); the Transformer, the SoccerMon case study and the simulations exist to argue that this table deserves that status. Classical exponentially weighted moving averages and related short- and long-horizon scalar summaries arise as restricted, typically univariate special cases of the same operator class. A simulation study with known memory mechanisms shows that \(Q\) and rotation-invariant recovery scores discriminate genuine multivariate or personalized memory from negative controls and from misspecified windows, whereas regime classification accuracy alone does not. SoccerMon serves as an empirical case study: a constructed Latent Memory Table attains \(Q\approx 0.73\) versus about \(0.40\) for classical and lagged principal-component baselines, with incremental held-out value for some wellness targets and Procrustes ensembles for row-wise reliability.

A sample of parametric Bayesian inference applications from astronomy, cosmology and particle physics is studied, augmented by mock data sets and toy problems. The parameter spaces of these parametric physical models and their posterior distributions from analysing specific data are characterized by (1) the number of model parameters, (2) whether the posterior shape is similar to a Gaussian, (3) whether the posterior has light or heavy tails, (4) how small the posterior is compared to the prior, i.e., how informative the data are, (5) whether some parameters remain unconstrained while others are highly constrained, (6) whether the posterior has multiple, disconnected modes, and (7) whether the inference undergoes phase transitions. These axis define a parameter space of inference problems. We characterize each of the inference problems and observe that inference in astrophysics spans the entire parameter space, from low to high dimensionality, mono- to multi-modal, and a variety of complex distributions that range from uninformative to highly informative. Furthermore, the computational cost of the physical models can range from milliseconds to dozens of seconds. The collated sample of inference problems is proposed as a standard test bed for new samplers. For reproducibility and ease of use, a Docker compute image is provided.

Bayesian model selection couples a discrete model indicator with a model-specific continuous parameter space. We introduce structured dimension-matched variational transdimensional inference (SM-VTI) for finite enumerable model spaces. A rooted construction graph expresses a model as a sequence of local stop/child decisions. Each typed edge compiles a declared scientific parent-child edit into an exact native-coordinate dimension-matching lifting; an edge-conditioned flow then learns the residual continuous transport. The resulting local policy and conditional flow define one direct joint variational distribution, without embedding every model in a saturated maximum-dimensional surrogate. We derive its exact path density and optimize the joint reverse-KL objective. On a controlled 15-model target, SM-VTI-Joint recovers terminal masses, local actions, and nonlinear conditional geometry. On a 128-model misspecified robust variable-selection problem, a 10-data-set nearly parameter-matched affine comparison with AVTI shows stronger early model-mass recovery and competitive final joint accuracy under the same target-evaluation budget.

Gradient-based Markov chain Monte Carlo methods are often introduced as a catalog of algorithms: Hamiltonian Monte Carlo (HMC), the Metropolis-adjusted Langevin algorithm (MALA), the No-U-Turn Sampler (NUTS), and several underdamped variants. This presentation obscures the common structure of the methods and, more importantly, the reasons why a sampler that is correct in principle may be ineffective in practice. We develop a unified account, beginning with exact continuous-time dynamics that represent idealized sampling methods and for which Metropolis adjustments are not required. Numerical discretization makes the dynamics computationally feasible but introduces bias. Metropolis adjustment removes the asymptotic bias by converting numerical errors into rejection, leading to HMC, MALA, NUTS, and the Metropolis-adjusted kinetic Langevin algorithm (MAKLA). The second half of the paper presents geometric design choices that determine practical performance, namely, although MAKLA and NUTS have nice theoretical properties, their sampling efficiency may be slow in practice. Importantly, a fixed mass matrix can whiten globally anisotropic targets, often fixing sampling inefficiency in Bayesian posteriors with large data. Whereas hierarchical posteriors introduce their own problem, causing state-dependent variation in the Hessian (e.g., Neal's funnel). We explain how a randomized step size can be used effectively to sample from such a distribution. The resulting paper is both a tutorial on the mechanics of gradient-based sampling and a set of practical recipes to improve sampler performance.

Metamodels are instrumental in reducing the computational burden associated with nested reliability analyses and optimization loops in Performance-Based Risk Optimization (PBRO) of structures under stochastic loads. In this context, stochastic emulators are particularly useful because they approximate response distributions while accounting for the intrinsic stochasticity of the simulator. Among these methods, Stochastic Polynomial Chaos Expansion (SPCE) is especially attractive because it does not require replications of nonlinear analyses at fixed input conditions. However, SPCE may present limitations in accurately representing extreme responses in the tails of structural response distributions. To address this limitation, this study proposes a framework that combines Generalized Stratified Sampling (GSS) with SPCE. The GSS scheme partitions the input space into strata according to the intensity of the hazard, improving the representation of extreme responses, while independent SPCE emulators are trained within each stratum. The conditional exceedance probabilities estimated in each stratum are then recombined using the total probability theorem to evaluate the probabilistic constraints. The proposed GSS-SPCE framework is applied to the optimal design of buckling-restrained brace cross-sectional areas in a two-story steel building. The objective is to minimize the initial construction cost while satisfying prescribed probabilistic performance constraints. Results show that the proposed framework accurately estimates structural response distributions, including their tail regions, while substantially reducing the number of nonlinear model evaluations required for PBRO.

Sequence-to-function models predict molecular phenotypes from complete sequence windows. At genome-wide association study loci, however, the prevailing design perturbs only the lead variant on the reference genome, even though the lead is often correlated with nearby variants through linkage disequilibrium. This single-variant perturbation implicitly fixes all linked alleles at their reference-genome states and may therefore create an uncommon or unobserved population haplotype. We study this input-construction problem at 38 Alzheimer's disease loci. We introduce HaploPerturb, which uses phased ROS/MAP genotypes or the European 1000 Genomes panel to fit a fixed-margin latent Gaussian factor model and rank partner configurations conditional on each lead allele. The leading public-panel construction agrees with the donor-panel construction at all loci under a strict linkage-disequilibrium threshold and 36 of 37 loci under a broader threshold after restricting to shared partners. Known-truth simulations show exact recovery of the dominant configuration under strong linkage disequilibrium and expose persistent residual correlation under a misspecified one-factor model. In an AlphaGenome benchmark against cell-type-specific ROS/MAP eQTLs, broader-set public-panel haplotypes yield microglial enrichment of 2.07 (95% whole-locus bootstrap percentile interval 1.48–3.60), compared with 1.43 (0.69–2.29) for a lead-only edit. Empirical-mode and LD-sign backgrounds yield 2.20 (1.60–3.64), with no detectable advantage or loss relative to the HaploPerturb top configuration. Thus population-informed sequence construction matters in this application, while the choice among reasonable leading haplotype rules is less consequential than the choice between a haplotype and a lead-only reference background.

We introduce a new class of uniformly ergodic MCMC algorithms, termed Diffeomorphic Contraction Sampler (DCS), and provide fast non-asymptotic mixing guarantees for DCS targeting distributions on $\R^d$ with arbitrarily heavy polynomial tails. DCS provides a solution to a well-known problem for MCMC samplers, which typically struggle with the combination of unbounded high-dimensional state space and vanishing gradients. The DCS pulls back a target on $\R^d$ onto a Euclidean ball $B(R)\subset\R^d$ and then samples from the transformed density on the convex set $B(R)$ via algorithms such as the Ball Walk, Hit-and-Run and others. A radial diffeomorphic contraction is chosen so that the pull-back density on $B(R)$ is bounded, implying uniform ergodicity for all targets with a finite polynomial moment. Non-asymptotic bounds for DCS require stronger assumptions such as log-concavity of the pull-back density. In practice, this is achieved approximately by a preconditioned automorphism of the ball $B(R)$, tuned via Variational Inference. Numerical simulation tests demonstrate that the DCS outperforms significantly the No-U-Turns sampler on multi-dimensional heavy-tailed targets arising as real-world posteriors in PosteriorDB benchmark. DCS also numerically outperforms in high-dimensional examples recently developed spherical projection samplers for heavy-tailed target distributions.

The effect size and nuisance parameters needed to appropriately size a clinical trial are generally not adequately understood before a trial begins. Through pre-planned early stopping rules, a conservatively planned group sequential design can effectively adapt the sample size for a clinical trial for a new treatment that is ineffective, very effective or minimally effective. One potential issue with such designs is that a substantial number of patients can be enrolled after a data cutoff for an interim analysis while data are being entered, cleaned, analyzed and discussed. A strategy of re-estimating the sample size at an interim analysis based on conditional power has been proposed to reduce somewhat this enrollment overrun issue. A purported advantage sometimes claimed is a smaller up-front planned sample size than a conservatively planned group sequential design. We demonstrate derivation of 2-stage group sequential and conditional power designs using the gsDesign R package and suggest comparing designs with comparable power using expected sample size calculations. While there are cases where conditional power designs may have small advantages, it is quite easy to derive very inefficient conditional power designs. This, along with the fact that a conditional power design may reveal something about the interim treatment effect, will often leave a group sequential design as the design of choice.

Scalable spatiotemporal modeling remains challenging because conventional methods rely on covariance models that tightly couple spatial representation and temporal inference, often leading to high computational costs. To address this difficulty, we developed coarse-to-fine spatiotemporal modeling (CF-STM), a framework that extends coarse-to-fine spatial modeling (CF-SM) to spatiotemporal settings. CF-STM represents latent spatial processes through multiscale locally weighted models, with temporal dependence modeled separately through state-space models defined at local centers. This covariance-free formulation achieves scalable computation while allowing temporal models to be flexibly specified without altering the spatial representation. Monte Carlo experiments demonstrate predictive performance comparable to that of alternative scalable space-time models at a substantially lower computational cost. An application to long-term residential land price data in the Tokyo metropolitan area shows that CF-STM flexibly captures complex spatiotemporal patterns while enabling interpretable inferences. CF-STM is implemented in an R package spCF (https://cran.r-project.org/web/packages/spCF/).

For a Markov kernel $T$ with an invariant probability measure $\pi$, we give a self-contained proof of the Markov chain convergence theorem via a criterion called asymptotic equivalence with the target. It assumes two parts about the Lebesgue decompositions of $T^n_x$ and $\pi$ for every starting point $x$: 1.) asymptotic absolute continuity: the singular mass $\mathrm{sing}(T^n_x \mid \pi)$ tends to $0$; and, 2.) asymptotic domination of the target: the singular mass $\mathrm{sing}(\pi \mid T^n_x)$ tends to $0$, as $n \to \infty$. Assuming a jointly measurable density for the absolutely continuous part of each iterate $T^n$ w.r.t. $\pi$, this criterion is sufficient and necessary for convergence. A positive minorant density version of it is verified in three cases: i.) $T$ has a positive transition density w.r.t. $\pi$; ii.) $T$ consists of an absolutely continuous part with positive transition density together with an atom at the starting point, which covers the Metropolis-Hastings algorithm; iii.) the transition density is positive only after a finite number of steps that may depend on the starting point $x$. To demonstrate our general criterion, we investigate the Gibbs sampler with random scan and the parallel tempering algorithm. Furthermore, we show that in all mentioned settings Birkhoff's ergodic theorem applies, so as to obtain the strong law of large numbers. Throughout this paper, neither irreducibility, nor aperiodicity, nor recurrence, nor couplings, nor splitting constructions, nor small sets are used. In all results, the state space is a general measurable space with no structure beyond a $\sigma$-algebra. That joint measurability is assumed of the Markov kernel, not of the space; countable generation supplies it. None of the theorems proved here is new; what is offered is a short route to a single, widely applicable Markov chain convergence criterion.

The conditional forecast-revision scale $\It=\{\Var(\E[X_{t+1}\mid\F_t]\mid\F_{t-1})\}^{1/2}$ measures the history-specific size of the forecast update induced by observing $X_t$. Because it is a conditional second moment built from two unknown conditional means, it is not directly observed. We study which estimator of $\It$ should be used under different structural assumptions and computational budgets. The comparison includes a block bootstrap, a conditional-variance model, a fitted state-space model, two $O(1)$ streaming smoothers, and the forget gate of an already-trained recurrent network. An error decomposition separates one-step-prediction error from conditional-second-moment tracking error. We show that externally tuned lag-only smoothers can be inconsistent when $\It$ changes at the sampling scale, although they attain the usual $T^{-2/3}$ mean-squared-error rate ($T^{-1/3}$ for $\It$) under slow variation; a correctly specified state-space estimator escapes this limit by using the current state. In volatility-driven designs, a cheap conditional-variance model is more accurate and over one hundred times cheaper as a point estimator than the implemented block bootstrap, whose value lies in the sampling distribution it provides rather than in point tracking. In state-driven designs, only the structurally matched filter recovers the fast variation. Read directly, a trained network's forget gate does not track $\It$ — though a supervised linear probe on the full gate vector does, so $\It$ is linearly decodable but not available for free. These results yield a practical rule: identify the conditional-second-moment structure, match the estimator to it, and then choose the least costly adequate method.

We develop a class of general diffusion-based stochastic particle optimisation methods for loss functions with intractable gradients. Specifically, we consider problems in which the loss gradient is an integral with respect to a parameter-dependent distribution, a structure that includes training generative models, fine-tuning, and learning latent-variable models. We introduce mean-field dynamics and its interacting-particle approximations, which contain several existing algorithms as special cases and provide a route to constructing new methods. Under well-posedness and joint contractivity assumptions, we prove exponential convergence and show that the continuous-time particle system admits a non-asymptotic error bound. We illustrate our framework by developing higher-order Langevin and momentum optimisers and evaluating them on maximum marginal-likelihood estimation and energy-based-model training.

We provide new estimates in Wasserstein distance for the asymptotic bias of the unadjusted Langevin algorithm, in the classical setting of log-smooth strongly log-concave measures. Our bound implies a Wasserstein mixing time of order $\kappa \sqrt{d}/\varepsilon$, where $\kappa$ is the condition number, $d$ is the dimension, and $\varepsilon$ is the target precision: this improves by a factor of $\sqrt{d}/\varepsilon$ over the previous state-of-the-art results.

Target trial emulation structures observational causal analyses around the protocol of an ideal randomized trial. By aligning eligibility, treatment assignment, time zero, and follow-up, it can reduce avoidable biases, but implementation still requires coordinated decisions about data construction, inverse probability weighting, diagnostics, outcome models, standardization, competing risks, and uncertainty estimation. This article provides a self-contained methodological guide and practical tutorial for TTE, an R package for target trial emulations with longitudinal observational data. We describe target trial protocols, intention-to-treat and per-protocol estimands, identification assumptions, baseline and person-period data structures, temporal ordering for longitudinal weights, stabilized treatment and censoring weights, weight truncation, balance and effective-sample-size diagnostics, weighted pooled discrete-time survival models, model-based standardization, competing-risk analysis, weighted Kaplan-Meier and Aalen-Johansen estimation, and cluster bootstrap at the original-individual level. Two fully synthetic examples illustrate end-to-end workflows: sodium-glucose cotransporter 2 inhibitor versus dipeptidyl peptidase-4 inhibitor initiation with all-cause death, and sequentially nested angiotensin receptor blocker versus calcium channel blocker trials with heart-failure hospitalization and competing death. The examples show how to obtain and diagnose estimates in R and interpret relative hazards, absolute risks, cumulative incidence, and differences between intention-to-treat and per-protocol effects.

We propose a robust deep mixture model based on a pathway-wise shared scale-mixture construction. Layer-specific component indicators are independently distributed according to their corresponding mixing proportions and jointly define a complete pathway through the latent hierarchy. Conditional on the selected pathway, a single gamma-distributed latent precision variable is shared across the deepest latent distribution, every intermediate latent transition, and the observation model. Integrating out this shared precision yields an exact multivariate Student-$t$ distribution for each complete pathway, allowing robustness to propagate coherently throughout the entire latent hierarchy rather than being introduced separately within individual latent layers. Model parameters are estimated using a stochastic expectation–maximisation algorithm. Complete-pathway responsibilities are evaluated analytically, whereas the shared latent precision variables and latent Gaussian variables are generated from their conditional distributions before updating the model parameters. The pathway-specific degrees-of-freedom parameters are estimated by one-dimensional numerical optimisation. Simulation studies demonstrate accurate recovery of the pathway-specific degrees-of-freedom parameters together with consistently improved clustering performance relative to the deep Gaussian mixture model under heavy-tailed and contaminated settings. Real-data applications further illustrate the ability of the proposed model to identify heterogeneous latent structures while reducing the influence of atypical observations. The proposed framework retains the hierarchical representation and parsimonious parameter-sharing structure of the deep Gaussian mixture model while providing coherent pathway-wise robustness.