Papers

Circular data recorded over time, such as wind direction, are increasingly observed in intervals rather than as single point-valued observations, yet no existing framework accommodates data that are simultaneously circular, interval-valued, and temporally dependent. Time series models for point-valued circular data, including recently proposed score-driven approaches, assume that each observation is a single angle. This paper addresses this gap by proposing a dynamic conditional score framework for interval-valued circular time series. The lower and upper boundary sequences of the observed intervals are treated as circular time series governed by score-driven dynamics, which preserves the periodicity of the circular domain and the dependency of boundaries, in which the boundaries are coupled through a toroidal distribution and a jointly derived score. Maximum likelihood estimation is discussed, and their finite-sample performance is evaluated through simulation. An application to interval-valued wind direction data illustrates the practical value of retaining the interval structure within a fully rotationally invariant circular framework.

Background: Credit risk modelling is fundamental in financial institutions, typically using credit history, repayment behaviour, income, and outstanding debt to estimate a borrower's probability of default. Although transactional data are sometimes incorporated through aggregated measures such as recency, frequency, and monetary indicators, behavioural information embedded in raw transaction narratives remains largely unexplored due to its unstructured nature. This limits opportunities to extract richer behavioural insights. Aim: This research investigates payment narratives from transactional data as an alternative data source for credit risk modelling. Setting: The study uses transactional data from a large South African financial institution focused on small and medium enterprises. Method: Techniques such as natural language processing and regular expressions are used to transform unstructured transaction descriptions into meaningful, explainable behavioural features by categorising transactions into interpretable spending categories. A novel categorisation algorithm was developed and evaluated for classification accuracy. The extracted features are incorporated into a credit risk scorecard to assess their predictive value. Results: The categorisation algorithm demonstrates high classification accuracy. Using these features in the scorecard resulted in a significant increase in Gini compared to a scorecard without them. Conclusion: Text-based transactional information can be systematically extracted, categorised, and operationalised to enhance credit risk model performance. Contribution: The study develops a novel, integrated, explainable methodology for transforming transactional text into behavioural features, introduces a structured RegEx-based categorisation approach, and demonstrates their predictive value within traditional scorecard frameworks.

Spatial prediction is a critically important task in data analyses that involve spatially or spatio-temporally-referenced data. It is commonly carried out using hierarchical models, which comprise a process model describing the latent process of interest, and a data model describing how the process is observed. The target quantity is the predictive distribution, which is the conditional distribution of the latent spatial process given the observed data. Unfortunately, the predictive distribution is rarely available in closed form and can be difficult to sample from. Here, we show how recent developments in guided diffusion modelling from machine learning can address a wide class of problems in spatial prediction when the likelihood function derived from the data model is known and differentiable, possibly after relaxation. To show the versatility of the approach, we apply it to both Gaussian and non-Gaussian process models and a variety of data models that describe direct partial observations, linear and areal aggregates, counts, right censoring, and area-exceedance functionals. Where possible, we validate the framework against exact and Markov-chain-Monte-Carlo (MCMC)-derived predictive distributions. We demonstrate our approach on several simulation studies and in a controlled experiment of methane-flux inversion over the UK and Ireland. We show that guided diffusion can sample from the correct predictive distribution in a variety of settings, that it can be applied to models where MCMC is intractable, and that it scales well with the dimensionality of the latent process.

Determinantal point processes (DPPs) model subset selection probability with repulsion, but many applications also exhibit attraction arising from pairwise relations or shared latent structures. Nonsymmetric DPPs (NDPPs) capture both effects, but their kernels are difficult to interpret and specify. We propose the determinant-ratio point process (DRPP), which models repulsion and attraction via a ratio of the determinants of two symmetric positive-definite kernels. This formulation preserves the feature-based interpretation of standard DPPs while allowing the two effects to be specified separately. Furthermore, we introduce an auxiliary-variable representation that enables efficient MCMC sampling and yields tractable expectation formulas for key statistics including the normalizing constant, which together facilitate maximum-likelihood learning. We theoretically characterize the interaction structures well suited to DRPPs: for clustered attraction with $C$ latent groups, DRPPs exactly realize the target pairwise interactions with an attraction rank $C$, whereas NDPPs can require a rank that grows with the number of items given the same repulsion component. Conversely, nontransitive pair-specific attraction may require compensation through the DRPP repulsion kernel. Experiments support these predictions, showing that DRPPs fit transitive attraction with substantially fewer parameters than NDPPs.

Bayesian simulation-based inference (SBI) can struggle to recover increasingly concentrated posteriors as observation size grows. We consider learning sufficient summary statistics jointly with the posterior as a general principle for scalable inference in various SBI methods. A common measure taken in SBI is dimension reduction through predetermined summaries, which may partially discard posterior information. Our approach instead imposes structure on how the posterior learner depends on a learned summary (e.g., affinity), and learns the summary simultaneously with the posterior. We identify an important mechanism named two-way guidance: under a structural restriction on the summary, the posterior sufficiency of the learned summary simplifies posterior learning, whereas a suitable posterior objective guides the learned summary towards posterior sufficiency. This mechanism explains a counterintuitive finding that the structural restriction on the summary statistics improves posterior learning, rather than limit its performance. We establish theoretical guarantees on both posterior estimation error and the learned summary's departure from sufficiency, thus identifying conditions on the simulation budget for recovering increasingly concentrated posteriors as observation size grows.

Aggregating independently computed sub-sample inferences is a well-known challenge in large-scale Bayesian computation, federated learning, and opinion pooling. The core problem is to approximate a distribution whose density is proportional to a product of individual densities. We propose Product Monte Carlo (PMC), an algorithm for computing arbitrary expectations under a product of densities using only Monte Carlo draws from individual densities. We turn the intractable product expectation into a regression problem by expressing it through a family of computable conditional expectations. Besides the point estimate, our fully Bayesian regression formulation provides coherent uncertainty quantification and principled diagnostics. We establish theoretical guarantees on consistency and characterize convergence rates for the PMC estimator. Across simulated and real-data experiments varying the number, dimensionality, and heterogeneity of the component distributions, PMC achieves higher estimation accuracy, better-calibrated uncertainty coverage, and improved scalability at a lower computational cost than existing methods. In addition, PMC requires no density evaluation, making it applicable to simulation-based inference, as we demonstrate in an application with protein images.

We detail scalable methods for approximating the quadratic entropy $p^T d p$ for arbitrary distributions $p$ and common distances $d$ of negative type. We focus on the Euclidean and spherical geodesic cases, which both use random feature embeddings and projections to dramatically improve computational complexity within a simple framework. Amortization of a single large matrix multiplication and control variates further enable computation at large scale with low memory and runtime in situations where $d$ is held constant while $p$ varies. We demonstrate this with a comparison against direct pair sampling and bibliometric/scientometric examples on Open Graph Benchmark datasets, revealing papers, fields, and institutions with both particularly narrow and broad interdisciplinary reach from their citations and text features alone.

We present the postshock R package, which implements and extends a donor-based framework for forecasting when a structural shock is known and the target response of interest is not yet observed. The package estimates shock effects from historical donor episodes, balances donors using specified matching features, and transfers the resulting adjustment to a target-series forecast. It provides integrated workflows for conditional mean forecasting through ARIMA and ARIMAX models and conditional variance forecasting through GARCH-X models. Additional functionality includes structured donor pools, control-shock regressors, automated GARCH-X order selection, processed data objects, and reproducible empirical workflows.

Estimating covariate-dependent intensity functions from spatial point pattern is challenging when the integrated intensity function cannot be evaluated analytically. This challenge becomes even worse in Bayesian inference due to repeated integration of the intensity function over the high-dimensional covariate space. This paper develops a Bayesian soft-tree Poisson process model that combines adaptive covariate partitioning with smooth intensity function. The key idea is to use soft gating functions to smooth transitions between partition boundaries in a tree generating process while preserving analytic evaluation of integrated intensity. A Bayesian CART prior is adopted to learn the partition structure from data. Given the partition tree, conjugate priors are assigned on leaf parameters to enable analytic marginal likelihood computation. This paper also develops a reversible-jump Markov chain Monte Carlo (RJ-MCMC) algorithm for posterior inference. Posterior consistency is established for the proposed model with respect to the Hellinger distance. Extensive simulation studies and an application to wildfire data confirm the advantages of the proposed method over existing methods.

We provide the first generic exact sampling algorithms with finite expected total cost for uniformly elliptic multivariate diffusions and the corresponding diffusion bridges with bounded coefficients and bounded continuous first derivatives. We formulate a computational cost model that counts work in terms of generation of Gaussian and uniform random variables to execute acceptance–rejection sampling, as well as function evaluations of the drift and diffusion coefficients, together with elementary operations such as sums and multiplications. Our constructions combine randomized parametrix expansions with a finite weighted graph that incorporates both bridge endpoints into the proposal. The exact sampling of the diffusion is minimax optimal jointly in the horizon and number of requested observations, while bridge sampling is minimax optimal in the number of observations for fixed horizon and endpoints.

Small area estimation (SAE) often requires both borrowing information across areas and benchmarking estimates to reliable aggregates. We develop a Bayesian framework that addresses these two objectives jointly through a new family of Benchmarking priors. The priors are induced by a benchmark-constrained regularization problem. The resulting family includes a Benchmarking Prior that incorporates the benchmarking restrictions without additional regularization across areas, and Single and Multi-View Laplacian Benchmarking Priors that introduce regularization through graph Laplacians constructed from area similarities based on external covariate information. For posterior computation under these degenerate priors, we develop tailored MCMC algorithms based on a reduced parameterization of the constraint space. We assess the proposed models using a data-based simulation and apply the framework to estimate Average Household Size (AHS) at the municipality level in Colombia in 2025. In this application, graph-based regularization improves model performance, with the Multi-View models generally producing more precise municipality-level estimates.

Continuous structural connectivity describes the connection between any two points of the cortical surface by an intensity function. Based on the endpoints of tractography streamlines, existing methods estimate this function mostly by kernel smoothing, which only gives a point estimate. We propose a Bayesian framework that estimates continuous connectivity with uncertainty quantification. The endpoints of tractography streamlines, mapped to two spheres, are modeled by a log-Gaussian Cox process on the product of the spheres with a Matérn Gaussian process prior. With the finite element representation of the prior, the Laplace approximation of the posterior enables the computation on a fine mesh with over 400,000 latent variables per hemisphere pair. Numerical experiments show the accuracy of the point estimate of this approach and the coverage of its credible intervals. In the Adolescent Brain Cognitive Development Study, the proposed method predicts held-out streamlines better than kernel smoothing.

In statistical models for categorical responses, including binary, nominal and ordinal responses, a phenomenon called “separation” may occur which can render maximum likelihood estimate nonexistent and/or nonunique. Based on Sablica et al. (2026)'s general characterization of separation phenomena via structure vectors, in this article we discuss computational, numerical and practical aspects of separation phenomena in different software packages, including consequences, remedies and ways of diagnosing them. We further introduce the R package divoRce which provides the most comprehensive software suite to diagnose separation in existence. It allows to check for separation, to distinguish between types of separation, characterize the recession cone and identify variables and observations contributing to separation for almost any categorical response model in the literature, including models with baseline-category link, cumulative link, adjacent-category link and sequential link specifications. Diagnostics are based on methods for computational geometry and linear programming and can be done in floating-point or rational arithmetic with different solvers. The functionality is easily extendible by wrappers and S3 methods to accommodate implementations in third-party packages or of new categorical response models. We discuss our implementation in detail and show our software in action in a wide range of applications. We recommend to incorporate our checks for separation and existence of the MLE as a diagnostic routine in every categorical data analysis both in third-party software implementations and as part of a standard model diagnostic workflow.

In the analysis of multivariate data, factor models represent a powerful technique for both reducing dimensionality and facilitating qualitative understanding of the latent determinants governing the observed data. Interpretability, however, is often hindered by the non-identifiability of factor loading matrices due to rotational invariance. A prevailing strategy to address this issue is to consider positive lower triangular structures. A relaxation of the latter condition has been recently proposed in the literature with the so called generalized lower triangular structure. However, its Bayesian implementation presents substantial challenges, requiring complex reversible-jump algorithms and making prior elicitation difficult because of the limited interpretability of the model parameters. In this paper, we introduce a generalized lower triangular process prior that allows both global and within-component sparsity structures for learning both the number of factors and possible sparsity structures within the vector of observations. The proposed approach provides straightforward prior parameters elicitation exploiting possibly different prior information on the rank and sparsity characteristics. We also explore connections with the Indian buffet process, providing further insight into the sparsity structure induced by the proposed prior. Posterior computation can be performed resorting to a Gibbs sampler with Metropolis-Hastings moves. The efficacy of the proposed method is demonstrated through comprehensive simulations and the analysis of the Big Five Personality Test data.

We study the problem of sampling from a target distribution $\pi(x)\propto e^{-f(x)}$ supported on a convex set $ X\subseteq\mathbb R^d$, when the potential $f$ is accessible only through a stochastic first-order oracle. Mirror Langevin algorithms provide a natural approach to constrained sampling by transporting the problem to an unconstrained dual space and discretizing the resulting Mirror Langevin diffusion. Existing implementations, however, typically use a fixed discretization step size and consequently retain a nonvanishing discretization bias at any fixed step size. Moreover, their direct extension to settings with noisy gradient information entails the challenge of controlling both discretization and stochastic-oracle error. We study a stochastic first-order version of the Mirror Langevin Algorithm (sFO-MLA) and, as our main contribution, develop a warm-started two-loop implementation in which an outer loop progressively decreases the step size while an inner loop runs sFO-MLA (with a fixed step size) for an appropriately chosen epoch length. The construction provides a principled schedule linking step sizes and epoch lengths, so that successive epochs warm-start from increasingly accurate distributions rather than repeatedly paying the cost of mixing from a cold start. We establish finite-time Wasserstein guarantees for sFO-MLA that explicitly separate mixing, Euler–Maruyama discretization, and stochastic-gradient errors. These bounds yield a fixed horizon rate of $\widetilde O(T^{-1/2})$ and show that the two-loop scheme removes the associated logarithmic penalty, attaining the canonical $O(T^{-1/2})$ rate under a geometric step-size schedule and corresponding epoch lengths. We illustrate the methodology in two statistically distinct settings.

To address the critical need for statistical methods that evaluate practical magnitude alongside statistical significance, we introduce a formal graphical inference framework that intrinsically combines effect size with uncertainty. Built upon resampling and global envelopes, this approach provides universality and direct visual interpretability. While existing envelope tests provide only weak family-wise error rate (FWER) or false discovery rate control, we advance the methodology by introducing novel step-down global envelopes. We theoretically prove that several of these envelopes achieve exact strong FWER control under exchangeability, guaranteeing rigorous inference for each individual local hypothesis. Beyond yielding adjusted $p$-values, the framework provides adjusted subset $p$-values for blocks of hypotheses, such as the functional effect of a covariate within a specific category. By explicitly visualizing the size and direction of the effect relative to its variability under the null hypothesis, our approach overcomes the dichotomous nature of traditional testing and provides profound informational value. The proposed framework is applied to multivariate and functional linear models.

The logarithmic multiplicative error model (log-vMEM) has been useful in modeling and forecasting multivariate positive-valued financial time series. The number of parameters grow rapidly with the dimension of the system and the lag order, making estimation computationally demanding in high-dimensional settings. This paper describes regularized estimation via hierarchical lag structures (Nicholson et al., 2020) for log-vMEM models with multivariate gamma error distribution of Tsionas (2004). The parameter estimation is performed using a blockwise coordinate descent algorithm with a Gauss-Seidel-style update scheme (Wright, 2015). This enables an efficient computation strategy compared to traditional penalized maximum likelihood approaches. The competing models are juxtaposed against each other by combining three hierarchical lag structures (componentwise, elementwise, own-other) and four penalties(group-lasso, adaptive group-lasso, group-mcp, and group-scad). Extensive simulation runs have been performed to test the parameter recovery for both the unpenalized and the penalized models. We apply the proposed methods to model the joint dynamics of robust intraday realized volatility measures for Microsoft (NASDAQ: MSFT) for the competing models. The numerical integration step of the log-likelihood is identified to be the principal computational bottleneck. We address this issue by using GPU-accelerated quadrature integration thus improving computational scalability of the proposed models.

Predicting the shapes of entire curves requires accounting for nonlinear geometry and unknown alignments between curves. We develop elastic kernel ridge regression, a nonparametric method for planar curve responses with scalar, multivariate, or functional covariates. Using square-root velocity representations, we formulate penalized conditional Fréchet mean estimation in a shape space invariant to translation, rotation, scaling, and reparametrization. A vector-valued reproducing kernel Hilbert space provides a flexible nonlinear link through the spherical exponential map. We use an alternating algorithm to align each observed curve to its current fitted value, and updates the regression function. To facilitate this, we provide a new Euclidean quasi-Newton solver that exploits scale invariance of the reparametrization objective; this accelerates alignment while retaining accuracy in a numerical comparison. Simulations demonstrate the benefits of estimating alignment within the regression and respecting spherical geometry. Applied to vocal tract contours from a real-time magnetic resonance imaging recording, the method recovers missing frames and reconstructs shape trajectories from downsampled data. A speech inversion proof of concept on the same recording predicts tongue shapes from acoustic features, illustrating the method's potential. The method is implemented in the R package sphereg2.

Finite mixtures of Clayton, Gumbel, Frank and Gaussian copulas can describe dependence that differs between the upper and lower tails. We study a two-stage Bayesian analysis in which ranks or kernel estimates replace the margins and the copula likelihood is evaluated at the resulting pseudo-observations. This pseudo-posterior is strongly consistent for the copula density whenever the log density admits a logarithmic boundary envelope. Both stages may use the same data, margins may be standardized within observed strata, and neither smoothness nor identifiability is required. The envelope holds for finite mixtures of Gaussian, Student, Clayton, Gumbel and Frank copulas in any fixed dimension; tail-dependence coefficients and conditional tail probabilities are therefore consistently estimated. We further prove that Clayton, Gumbel, Frank and Gaussian copulas are jointly finitely linearly independent, which makes mixture weights and components identifiable and consistently estimated. In simulations the pseudo-posterior matches multi-start maximum pseudo-likelihood in large samples. It is more stable in small samples and near independence (every component close to the independence copula), where tail coefficients are learned long before the weights and a marginal Metropolis sampler is up to twice as efficient as data augmentation. Two physical fitness datasets show mirror-image asymmetries. Among 8772 university students, sprint and jump performance are coupled mainly at the top. Among 5336 adults in a national health survey, low grip strength and low daily activity cluster together, increasingly with age, whereas high values do not. A Gaussian copula misses both patterns in held-out data.

The Multilevel Markov chain Monte Carlo algorithm accelerates Bayesian inversion by exploiting correlations between discretized quantities of interest across a hierarchy of model resolutions. A key practical parameter is the subsampling rate used to generate coarse-level proposals for the finer-level chains. Existing practice sets this rate equal to an integrated autocorrelation time as a proxy for proposal decorrelation, but this choice can be overly conservative and lead to substantial computational overhead. In this work we treat finite subsampling as an error-budget question. We derive a cost-error criterion for selecting a computationally efficient subsampling rate using estimates of sample variances, integrated autocorrelation times, and per-level costs. The criterion balances multilevel sampling errors across levels rather than enforcing approximate independence of the coarse proposals. We clarify the theoretical role of independence through an ideal reset-reference kernel, for which independent coarse proposals may be replaced by proposals generated from a reversible coarse-level kernel. We derive a Poisson-equation representation of the perturbation bias resulting from finite subsampling. Under a state-averaged lower-level forgetting assumption, this yields a conditional geometric perturbation-bias bound; in computation, the corresponding projected diagnostic is used as an empirical error-budget check. Numerical experiments on a highly autocorrelated elasticity beam benchmark show substantial reductions in computational cost compared with IAT-based subsampling. A second Darcy-flow inverse problem gives supporting evidence of qualitatively similar finite-subsampling behaviour in another elliptic benchmark.

Spectral methods are central to high-dimensional statistics and machine learning, underlying procedures for covariance estimation, matrix denoising, representation learning, clustering, and latent variable modeling. In this work, we focus primarily on spectral estimators for high-dimensional factor models, where leading singular vectors of the data matrix are used to estimate latent structure and covariance parameters. In many modern applications, however, data are collected across related but heterogeneous tasks, studies, domains, or populations. Applying such estimators separately to each task can lead to unstable estimates when sample sizes are limited relative to dimension, while completely pooling the data can obscure meaningful task-specific structure. We introduce Hierarchical Spectral Shrinkage (HSS), a scalable empirical-Bayes framework for partially pooling task-specific spectral estimators. The method regularizes the leading spectral directions of each task toward a data-adaptivelylearned common basis, while allowing the amount of shrinkage to vary across tasks and spectral components. The proposed estimator arises as the posterior mean in a surrogate Bayesian regression formulation of the empirical singular vectors. For factor models, HSS yields partially pooled estimators of the group-specific loading spaces and covariance matrices and can be combined with different spectral estimators by replacing their empirical singular vectors with hierarchically regularized counterparts. More broadly, the same construction provides a mechanism for extending spectral estimators to collections of related but heterogeneous datasets. We demonstrate substantial improvements in estimation accuracy and out-of-sample performance over separate, fully pooled, and shared-subspace approaches in synthetic experiments and multi-study gene-expression data.

The log-rank test is the standard tool for two-sample survival comparison and has good power against proportional-hazards alternatives, but it loses power against sparse hazard departures, in which the hazard difference is concentrated in a small number of time intervals whose locations along the follow-up are unknown a priori (Kipnis, Galili and Yakhini, Biometrika 2026). Weighted log-rank tests – Gehan-Wilcoxon, Tarone-Ware, Peto-Prentice, Fleming-Harrington – each impose a pre-specified temporal emphasis. Combination procedures such as MaxCombo and the Yang-Prentice short-term/long-term hazard-ratio model relax that choice, but still scan only a small dictionary of global temporal shapes and remain insensitive to sparse departures occurring outside them. We introduce lifelines-hc, a Python package extending the lifelines survival library with the HCHG test: Higher Criticism applied to per-interval hypergeometric p-values. HCHG applies to right-censored two-sample data and detects sparse hazard departures at unknown locations without committing to any temporal pattern. Across four clinical case studies – CheckMate 057 PFS (n=582), COMET-1 OS (n=1028), AZURE DFS (n=3359), and the Copenhagen Study Group for Liver Diseases cirrhosis trial (n=446, fully public individual patient data) – HCHG achieves p <= 0.014 in every dataset, while the log-rank test is non-significant throughout (p >= 0.26). MaxCombo and the Yang-Prentice adaptive log-rank test detect the delayed-benefit crossing in CheckMate 057 (p < 0.001) but are non-significant on the remaining three (p >= 0.24), where the departure is sparse or multi-window rather than a smooth short-term/long-term hazard-ratio pattern. lifelines-hc is freely available under the MIT license at https://github.com/alonkipnis/lifelines-hc and via pip install lifelines-hc.

Statistical inference in sequential mediation models requires evaluating the cumulative distribution function (CDF) of the product of three normal coefficient estimators, an operation that entails integrating over a non-convex region with boundary $xyz=v$. While the Delta method provides an instantaneous first-order approximation, its asymptotic variance collapses whenever two or more path coefficients approach zero, producing severe undercoverage near parameter boundaries. Nonparametric bootstrapping avoids gradient collapse but incurs an $O(B \times N)$ computational cost that becomes burdensome in large-scale simulation studies or iterative power analyses. We propose a model-based dimension-reduction algorithm that integrates out the third variable analytically under the trivariate Gaussian distribution for arbitrary mean vectors and positive-definite covariance matrices. Partitioning the $xy$-plane into quadrants isolates the sign change at the coordinate axes, reducing the problem to an adaptive two-dimensional quadrature whose fixed-grid cost is $O(n^2)$ in place of $O(n^3)$. In simulation benchmarks across six parameter regimes, the algorithm achieves a mean absolute error of $3.1 \times 10^{-5}$ relative to a $10^{8}$-sample Monte Carlo reference, evaluating the distribution in under one second and yielding plug-in quantile confidence intervals whose empirical coverage was near or above nominal.

Sensing technologies have advanced rapidly across industries ranging from energy to automotive manufacturing. These systems generate high-dimensional (HD) data characterized by complex nonlinear patterns and strong temporal dependencies. Traditional statistical monitoring methods are often limited in their ability to capture such nonlinear structure. Likewise, many analytical approaches used in End-of-Line testing rely on predefined thresholds and heuristic rules, which restrict their ability to detect informative anomaly signatures in HD temporal data. In contrast, while modern deep learning and generative AI models offer strong predictive capabilities, they are often unsuitable in applications where data are costly to collect and where the monitoring system must remain interpretable, low-latency, computationally efficient, and usable by non-technical practitioners. To overcome these limitations, we propose an advanced multivariate monitoring framework for HD data. The framework operates in two stages. In the first stage, the data are preprocessed to remove incomplete and non-informative samples and to temporally align time series data. In the second stage, nonlinear dimensionality reduction is performed, followed by anomaly detection through a control chart based phase I monitoring procedure. The framework can be used in both unsupervised and supervised settings, depending on the availability of ground truth labels during training. Moreover, its flexible and modular structure allows practitioners to adapt its components to different domains and operational requirements. We evaluate the proposed framework on real production data from an automotive manufacturing environment at Ford Motor Company. The proposed method achieves higher accuracy, recall, and F1 score than the company's existing model, improving these metrics from 0.50, 0.30, and 0.429 to 0.625, 1.00, and 0.769, respectively.

In the context of functional data analysis, we propose a Bayesian test to assess whether the mean function of a Gaussian population is identically zero. By leveraging the structure of the reproducing kernel Hilbert space (RKHS) associated with the covariance function of the underlying process, we approximate the mean by a finite linear combination of kernel sections. This finite-dimensional representation arises from the natural assumption that the mean lies in the corresponding RKHS, enabling a parametric Bayesian approach to hypothesis testing through a spike-and-slab prior on the coefficients. Given the absence of a canonical reference measure on infinite-dimensional spaces, the likelihood is expressed via a Radon-Nikodym derivative between induced Gaussian measures. Since the number of kernel sections is unknown, we employ reversible jump Markov chain Monte Carlo (RJMCMC) to explore the trans-dimensional parameter space and estimate the posterior distribution. From the resulting posterior samples we construct estimators of the posterior probability of the null hypothesis and the Bayes factor, providing quantifiable evidence against the null and supporting informed decisions under uncertainty. Furthermore, to validate this methodology we establish near-optimal posterior contraction rates for the proposed test under mild conditions, demonstrate its empirical performance on simulated data sets, and illustrate its application with real climate data.

Multivariate functional data arise in many modern manufacturing systems, where multiple sensors record densely sampled process trajectories. Monitoring such data is challenging because nominal variation is strongly correlated across samples, sensors, and time, while faults may appear either as isolated deviations or as structured departures concentrated within a limited number of sensor-specific temporal trajectories. We propose two unsupervised sparse tensor decomposition methods that preserve this multimode structure. Entrywise Sparse CP Decomposition (ES-CP) uses an entrywise \(\ell_1\) penalty to identify localized anomalies, whereas Fiberwise Sparse-Group Lasso CP Decomposition (FG-Lasso) combines entrywise and fiberwise penalties to detect both localized deviations and anomalies concentrated within temporal fibers. Both methods represent nominal process behavior through a low-rank CP decomposition and are estimated using alternating optimization with closed-form sparse-component updates. Two simulation studies evaluate performance under different fault structures, signal severities, noise levels, and missing observations. FG-Lasso attains or ties the highest macro F$_1$ score in almost all settings in the first study and achieves the highest macro F$_1$ score. In a multichannel forging-process case study, FG-Lasso and ES-CP obtain macro F$_1$ scores of 0.85 and 0.82, respectively, compared with 0.69 or lower for TRPCA and PCA-based anomaly detectors. The results demonstrate that explicitly matching the sparse penalty to the anticipated fault structure improves both anomaly detection and fault localization in high-dimensional functional processes.

We introduce an algorithm for efficiently preconditioning log-concave and log-smooth distributions that scales logarithmically with the condition number of the underlying distribution. Based on Gaussian cooling, our multistage method approximately samples from a sequence of well-conditioned distributions to construct a preconditioner for each subsequent stage of the cooling schedule. This method is motivated by existing practical approaches for mass matrix estimation in Markov chain Monte Carlo (MCMC), in which an effective preconditioner for the target distribution is estimated from samples produced during a warm-up period. However, our approach is conceptually distinct from existing approaches and moreover permits theoretical analysis. Our algorithmic framework is flexible, allowing essentially any log-concave sampling algorithm to act as the subroutine within each cooling stage. For instance, using first-order rejection sampling (FORS) as the sampler, the first-order query complexity of our method is $\tilde{O}(d^{4/3} \log \kappa)$ for a large class of log-concave and log-smooth distributions with condition number $\kappa$ and $\tilde{O}(d^{1/3} \log \kappa)$ for translation-invariant distributions. In the translation-invariant setting, we provide numerical experiments for the lattice $\phi^4$ model, a simple lattice field theory.

Sampling from a probability distribution with access only to its log-density is a foundational problem in statistics, machine learning, and the physical sciences. This problem is studied in the data-free score-based diffusion literature, which samples the target density via estimating scores of densities produced by convolution with Gaussian random variables and implementing a reverse SDE. In this framework, score estimators are constructed via Markov chain Monte Carlo, importance sampling, or rejection sampling. The estimators often operate far beyond the regime where mixing time and error control can be established and have large bias and variance when estimating the underlying score. In our method we propose post-processing: splitting a fixed computation budget across multiple evaluations, querying the estimator at different points, and fitting a single pooled estimate via linear combination to minimize bias and variance at a target location. The final value is an affine map of a vector of multiple score estimates. We show this reduces score estimation error when compared to the computation-matched base estimator queried at a single point. We demonstrate improved sample quality over other score estimation methods matched for computation cost on difficult target densities such as log-concave densities with large condition number, multimodal densities, and high-dimensional densities.

The log predictive density score (LPDS) and related estimators of the expected log predictive density (ELPD) have become the de facto standard for Bayesian model comparison. The LPDS is a natural out-of-sample measure of predictive density performance that accounts for parameter uncertainty, is less dependent on the prior than the marginal likelihood, and is easily computed sequentially using standard Monte Carlo from posterior draws. Our article demonstrates, however, that the standard Monte Carlo estimate of the LPDS can be extremely variable and severely biased when the test data has outliers, leading to misleading conclusions and suboptimal model choices. Estimating the LPDS as the difference between the log marginal likelihood of the full dataset and the training data is shown to be a dramatically more efficient approach. This estimator is also shown to be a much more effective fallback for observations that have been flagged as problematic by the widely used PSIS-LOO model comparison framework in Stan.

We consider the computation of the network modularity, a common measure of the structure of networks or graphs that quantifies the division of a network into modules or clusters. The computation of the network modularity is an NP-hard problem, thus making its exact computation often infeasible in practice. This has resulted in the development of several heuristics in the literature. In this work, we consider the exact computation of the network modularity with the help of a branch-and-bound algorithm. Our algorithm is guaranteed to find the optimal clustering that maximizes the network modularity, however it achieves this without a full exploration of the search space. This is accomplished by bounding partial solutions and discarding them if it can be foreseen that a partial solution will not yield an improvement over an existing solution. We assess our algorithm with respect to accuracy and runtime and compare it to a brute-force approach as well as two state-of-the-art benchmarks, simulated annealing and the Louvain method.