Papers

When data are grouped, hierarchical or multilevel models are commonly used to account for group-level variation with group-specific parameters. Leave-one-group-out cross-validation (LOGO-CV) is a suitable tool for evaluating predictive performance for new groups, providing an estimator of the expected log predictive density (elpd). Brute-force LOGO-CV requires one model refit per held-out group, often using computationally expensive inference algorithms such as MCMC. This is costly, particularly for large numbers of groups or complex model structures. Commonly used importance sampling approximations, intended to reduce this cost, tend to fail because the group-specific parameters of the held-out group must be integrated out. We identify two key challenges in LOGO-CV elpd estimation: approximating the LOGO posterior and computing the grouped marginal likelihood. We compare 11 strategies, including 5 newly proposed, to address them. Among others, we combine Pareto-smoothed importance sampling or adaptive importance sampling with integration techniques such as Laplace approximation, adaptive Gauss-Hermite quadrature, and bridge sampling. We evaluate these strategies in both simulation experiments and real-world case studies, which show that marginalising over the group-specific parameters substantially improves the reliability of the importance sampling approaches.

Single-cell RNA sequencing studies routinely use clustering to define putative cell types and cell states, yet the observed separation may arise from sampling variability rather than genuine biological heterogeneity. This paper studies formal significance testing of such clustering structure in high-dimensional data. Existing SigClust methods assess clustering significance through Monte Carlo simulation under a Gaussian single-cluster null, but this assumption can be unreliable for normalized gene expression data and other non-Gaussian settings. We propose SigClust-LCP, a nonparametric extension that models a single cluster by a log-concave distribution. To make this approach computationally feasible in moderate to high dimensions, we develop a score-matching estimator for log-concave projection inspired by recent generative modeling ideas. We establish theoretical guarantees for the estimator and for its use in clustering significance testing. Simulations show that SigClust-LCP controls Type-I error more reliably than existing methods across a range of unimodal and mixture distributions while retaining competitive power. In a single-cell RNA sequencing analysis of Hydra cells, the method avoids spurious subclusters within annotated cell populations and supports biologically meaningful separation across lineages and body-axis regions.

We develop a finite-horizon calibration method for raw-score martingale posteriors, with von Mises–Fisher models as the main worked example. Starting from the maximum likelihood estimator, predictive paths are generated by simulating future observations from the current fitted model and updating the natural parameter by unpreconditioned score increments. The main methodological step is to separate predictive simulation from covariance calibration. Raw-score increments have Fisher-information covariance, whereas Bernstein–von Mises calibration requires inverse-information covariance. We therefore apply a terminal linear correction based on a local information estimate. For more efficient implementation, we also introduce a hybrid version that replaces the omitted tail of the infinite predictive continuation by a Gaussian approximation with matching leading-order quadratic variation. We prove fixed-n convergence, a finite-horizon approximation bound for the Gaussian tail, and a Bernstein–von Mises limit for the hybrid post-corrected sampler under local regularity and consistent terminal calibration. Simulations show that tail correction reduces truncation-induced underdispersion, and an OSCAR ocean-current example illustrates local directional uncertainty summaries.

Cross-validation is a resampling procedure that provides a point estimate of generalization error for any predictive model. Cross-validation is widely used for model selection and evaluation. Uncertainty in the cross-validation estimate is challenging to quantify, and estimation of its variance is known to require multiple runs of the entire resampling procedure. Nested cross-validation computes a prediction interval for the generalization error for a given model and training dataset by resampling the entire cross-validation procedure but incurs extraordinary computational cost. We provide an efficient method for computing the nested cross-validation prediction interval using only a single model fit for some penalized regression models including ridge regression, spline smoothing, and some functional regression models. We characterize when our proposed method should be expected to out-perform resampling-based nested cross-validation in various scaling regimes as well as in finite samples. Experiments for functional principle components regression demonstrate non-trivial cases in which our proposed method improves run times substantially.

We present a framework for the analysis and unmixing of random density mixtures in the Bayes Hilbert space. General identifiability results for mixtures in Hilbert spaces are established and applied to the Bayes Hilbert space setting. Building on these results, we propose a penalised maximum likelihood approach for the unmixing of Bayes Hilbert mixtures aimed at recovering the statistically space-efficient representation, together with a computationally efficient coordinate-wise maximisation algorithm for its implementation. The methodology is illustrated through a hyperspectral data application, where observations can be naturally embedded in the Bayes Hilbert space and analyzed in terms of distributional shape rather than amplitude. A complementary simulation study demonstrates the interpretability and practical performance of the proposed approach.

Process-based forest models are widely used to simulate forest growth and responses to environmental change, but their calibration and application often require many computationally expensive model evaluations. We present an implementation of the Physiological Processes Predicting Growth (3-PG) model in JAX that uses just-in-time compilation, vectorization, and GPU acceleration to reduce execution time. The implementation also supports automatic differentiation, providing gradients of model outputs and calibration objectives with respect to model parameters. This enables efficient gradient-based optimization and gradient-informed Bayesian calibration, extending 3-PG beyond conventional gradient-free approaches. The implementation produced results numerically consistent with r3PG for the evaluated configuration. Overall, the JAX implementation provides a faster and differentiable framework for calibrating and applying the 3-PG model.

Stochastic gradient Markov chain Monte Carlo (SGMCMC) methods enable scalable Bayesian inference, but their performance depends strongly on hyperparameters such as the step size, mini-batch size, and number of leapfrog steps. Since most SGMCMC algorithms lack a Metropolis-Hastings acceptance rate, standard acceptance-based tuning methods are not directly applicable. We propose HyperMC, a multi-fidelity tuning framework that combines Hyperband-style resource allocation with kernel Stein discrepancy (KSD) evaluation. By running multiple successive-halving brackets, HyperMC balances broad exploration of a continuous hyperparameter space with increasingly accurate evaluation of promising configurations under a fixed computational budget. We further introduce Robust HyperMC, which uses global grid initialization followed by elite-guided local refinement to reduce sensitivity to random candidate generation and noisy finite-budget evaluations. Under suitable approximation and concentration conditions for the estimated KSD, we establish that the successive-halving component selects a near-optimal configuration among the sampled candidates with high probability and derive a sufficient computational budget for successful selection. Experiments on logistic regression, probabilistic matrix factorization, and Bayesian neural networks show that HyperMC improves posterior approximation or predictive calibration relative to MAMBA, grid search, and heuristic baselines, while Robust HyperMC yields more stable and reproducible tuning results.

Understanding the generative mechanism of real-world networks is crucial for analyzing connection patterns and making inference from network data. Graphons are widely used to model such mechanisms. Network histogram methods are nonparametric approaches based on blockmodel approximations that provide an intuitive view of network connection structures. However, there is a lack of software packages that construct network histograms. We introduce the R package nethist, which implements network histogram-type graphon estimators within a unified interface. This package is applicable to both single-layer and multilayer networks, and includes graphical summaries for examining both local and global network structures. By providing comprehensive network analysis tools, nethist facilitates understanding of complex systems of interrelated vertices.

Multinomial probit choice probabilities over n alternatives are Gaussian orthant integrals, computed by simulation for thirty years, one expensive integral per alternative. Inversion, which is to say determining item attractiveness consistent with a prescribed choice probability vector, is even more difficult and has been considered impractical for correlated contests when n is large. Yet here, for families lying within a grammar including factor, block and hierarchical covariance structures, we exhibit a calibration tested at n = 1,000,000 reproducing probabilities to very high accuracy, even in the extreme tail. We must return to much smaller problems for any performance comparison to be possible due to limitations of the prior art. The Geweke-Hajivassiliou-Keane simulator is the standard (and still appropriate for high rank) but is two hundred times slower already at n = 200, and its measured cost grows roughly as n^2.8 while ours is linear. Furthermore our approach applies to any continuous performance distributions within reason: the Thurstone-Mosteller model families thereby become a practical alternative to logit at modern scale.

Sampling from distributions conditioned on desired semantic properties is an emerging challenge in modern generative modeling. Metropolis-Hastings (MH) provides a principled route to conditional sampling, but requires access to exact pointwise target-density evaluations, which are not available in generative settings. Meanwhile, pairwise comparisons by humans or model "judge" are highly accessible and have proved valuable across diverse applications. We introduce Pref-MH, a general exact MH sampler for judge-induced conditional distributions using only stochastic binary pairwise comparisons. Our key observation is that the MH unnormalized density ratio matches the preference odds of the Bradley-Terry (BT) choice model. The central challenge is that while MH requires precise ratio computation, BT judges provide only sampled binary feedback. To this end, we develop a valid accept/reject rule whose resulting Markov chain provably converges to the target distribution. We further show that, for a fixed proposal kernel and budget, Pref-MH is optimal in the Peskun-Tierney sense among this class of exact reversible acceptance rules. Experiments on text generation and molecular design with LLM judges, as well as image generation with VLM judges, demonstrate that Pref-MH provides a practical and flexible approach to conditional sampling when comparative feedback is relatively easy to obtain.

High-dimensional clustering is challenging when component distributions are both heavy-tailed and directionally asymmetric. We propose a deep skew-$t$ mixture model (DStMM), a hierarchical factor-analytic mixture based on the generalised-hyperbolic skew-$t$ normal mean–variance representation. A shared inverse-gamma mixing variable is propagated along each complete latent pathway, allowing heavy tails and directional asymmetry to be modelled jointly while preserving conditional Gaussianity. Each complete pathway therefore admits an exact GHST marginal representation. We formalise the reductions to symmetric deep $t$, Gaussian deep-mixture, and single-layer GHST factor-analytic models, discuss local non-identifiability and the implementation-level parameter-counting convention, and derive the conditional generalised inverse Gaussian law used for estimation. Estimation is carried out by a stochastic/Monte Carlo EM algorithm, with an explicit implementation-based parameter count for BIC architecture comparison. Simulation studies show that DStMM performs similarly to the symmetric robust model when skewness is absent but provides increasing gains as directional asymmetry becomes stronger, particularly under heavier tails; the same qualitative behaviour persists under smaller samples and unequal mixture proportions. Two real-data applications provide complementary evidence. On the UCI handwritten-digit benchmark, DStMM gives the strongest clustering performance under a common deep architecture, while on the Gas Sensor Array Drift data, DStMM improves on both deep Gaussian and deep $t$ alternatives and, under the implemented BIC criterion, selects a non-trivial second mixture layer. Together, these results support the value of propagating skewness and heavy-tail variation through a deep latent mixture while retaining an exact pathway-level likelihood.

Matrix-variate data with missing entries arise frequently in applications where observations are naturally organized as two-dimensional arrays. Although the matrix normal distribution provides a parsimonious model through its Kronecker covariance structure, standard EM estimation can be computationally expensive because arbitrary missingness patterns typically destroy this separability in the E-step. In this paper, we propose an efficient partial EM algorithm for matrix-variate normal data with missing entries. The proposed method updates the conditional mean and covariance of the missing component through coordinate-wise approximations, avoiding repeated inversion of pattern-specific covariance matrices and avoiding construction of the full vectorized covariance matrix. We further develop a specialized update for submatrix missingness, where the missing-block precision retains a Kronecker product structure, and the covariance update can be carried out independently in the row and column directions. Simulation studies show that the proposed methods substantially reduce computation time compared with exact EM while preserving nearly identical observed-data likelihood across a range of dimensions and missing proportions. A real-data application to hyperspectral image patches demonstrates that the proposed imputation strategy can be embedded within a matrix-variate mixture model for simultaneous imputation and clustering.

Gibbs sampling and coordinate ascent variational inference (CAVI) are two basic coordinate-wise methods for statistical computation. Recent analyses under strong log-concavity establish convergence rates for versions of these algorithms that update one uniformly selected block at each step. We extend both results to arbitrary fixed, strictly positive selection probabilities. The rates are governed by a selection-adapted convexity constant $\lambda^{\star}_{\theta}$, defined using the block-smoothness constants and the selection probabilities $\theta$. The same constant yields a contraction of relative entropy for the Gibbs sampler and a contraction of the mean-field objective gap for random-scan CAVI. The new bounds recover the uniform-scan results, and are never weaker than the naive comparison based on the smallest selection probability. They provide a principled way to adapt the scan to heterogeneous block geometry using curvature information.

We propose a generalized Bayesian clustering with regression model for unaligned longitudinal binary outcomes, motivated by seizure diary data from the Human Epilepsy Project. Seizure diaries are sparse, irregularly observed, and vary enormously across patients. A single fully-specified generative model tends to be either misspecified or computationally inefficient. We address the challenge by two strategies. We set up a regression by way of clustering as model-based clustering using a mixture model. For the latter, we take a generalized Bayesian perspective which replaces the full likelihood with a loss-based update using a generalized likelihood. We combine a trajectory similarity loss and a regression loss, so that clustering is informed by both trajectory similarity and the prediction of outcomes The trajectory similarity loss is constructed by representing each trajectory as an (empirical) distribution of subsequences, called reads, and then is defined based on the sliced Wasserstein distance between these empirical distributions. This loss allows alignment-free comparison of sequences that are irregularly observed or temporally misaligned, and it scales quasi-linearly in trajectory length. The regression loss is the negative log-likelihood of a probit regression. A prior on the cluster-specific parameters is defined by way of a Dirichlet process prior on the mixing measure.

Constructing confidence regions for stochastic gradient descent (SGD) ideally requires estimating the asymptotic covariance matrix, a severe computational bottleneck in high dimensions. Traditional cancellation-based batch means methods bypass this estimation but require inverting a sample batch covariance matrix. This introduces strict mathematical degeneracy when the parameter dimension exceeds the number of batches. To address this problem, we utilize equal batch size batch means method and propose a simultaneous, marginal-friendly framework. The proposed marginal statistics has a asymptotic Student's $t$-distribution, and eliminates the matrix inversion step, entirely circumventing high-dimensional degeneracy. To achieve valid simultaneous coverage, we present an algorithm utilizing wild bootstrap samples drawn from a statistic as a function of only the diagonals of the variance-covariance estimator, and to further incorporate the contribution of cross-dependencies, we introduce an efficient Quasi-Monte Carlo procedure utilizing a $t$-copula approximation. Additionally, we integrate a Lugsail variance estimator to aggressively correct finite-sample bias and under-coverage. The proposed methodology delivers interpretable, simultaneous hyper-rectangular confidence regions that are statistically robust, memory-efficient, and strictly scalable for high-dimensional inference. The theoretical results are supported by extensive numerical simulation analysis through various aspects of dimension, number of batches and error structure.

Markov chain pool decoding (MCPD) devised by Knill et al. (1996) identifies likely positive clones from noisy pooled-test results. The standard MCPD estimates clone-wise posterior probabilities using Gibbs sampling, but it may allocate excessive computational effort to low-scoring clones. This paper focuses on parallelizing MCPD on GPU architectures. Whereas the standard MCPD employs systematic-scan updates, we propose a score-weighted update scheme that updates high-scoring clones more frequently. We prove that the stationary distribution of the proposed Markov chain coincides with the target posterior distribution. To enable efficient GPU parallelization, we further incorporate the unbiased MCMC framework of Jacob et al. (2020) and employ a slot-refilling technique based on the arguments by Glynn and Heidelberger (1991) about the coupling of Markov chains. Experiments involving 1,298 clones, 97 pools, and three true positives demonstrate improved recovery compared with uniform decoders, while maintaining high overlap under high-noise conditions.

We develop a multifidelity method for fusing low-resolution simulations with computationally expensive high-resolution simulations, which are run infrequently and are therefore prone to bias. We formulate this fusion as a constrained optimization under missing-not-at-random (MNAR) selection bias. This formulation searches for the exponentially tilted high-resolution distribution that minimizes KL divergence from the biased baseline, subject to moment constraints derived from low-resolution simulations. This optimization requires first estimating the biased baseline conditional density $f$ as a nuisance parameter. We estimate $f$ using a score-based diffusion model. To eliminate the generative model's regularization bias that harms the downstream task, we apply targeted maximum likelihood estimation (TMLE). TMLE debiases $\hat{f}$ via a targeted exponential tilting, rendering the target parameters insensitive to first-order nuisance estimation errors. To execute this computationally, we adapt generative model steering, a technique originally developed for human-preference alignment. Using Feynman-Kac steering with a reward function based on our formulation, we simultaneously execute the exponential tilts for MNAR and TMLE at inference time, avoiding expensive retraining costs. Code available [here](https://github.com/Jong-Min-Moon/multifidel_emul_by_FK).

We introduce a robust clustering method, MK-means DPD, that estimates cluster centers and covariance matrices using density power divergence (DPD) measures combined with Mahalanobis distance, making it resistant to outliers and adaptable to heterogeneous, elliptical clusters, unlike the classical K-means algorithm. Since Mahalanobis distance-based K-means lacks a general convergence guarantee, we further introduce a convergent variant, Density-Consistent MK-means DPD (DC-MK-means DPD), which redefines the cluster assignment step in terms of a pointwise DPD loss. We prove a formal theorem establishing that the resulting algorithm converges in a finite number of steps. We also propose two new robust internal evaluation indices, a Median Davies-Bouldin Index and a Trimmed Calinski-Harabasz Index, to ensure that performance comparisons are not themselves distorted by outliers. The efficacy of the proposed methods is demonstrated on simulated data, showing superiority over existing methods, and on two real datasets: Iris data, to identify similar species, and COVID-19 case fatality rate and infection rate data for countries worldwide, examining the resulting clusters' geographic and socio-economic patterns.

Kernel density estimation (KDE) is one of the most fundamental statistical estimators of density functions. Its direct implementation on a dataset of $N$ points incurs an $\mathcal{O}(N^{2})$ computational cost, which is prohibitive for large-scale datasets. Kernel approximation techniques can be applied to bring the computational cost down to $\mathcal{O}(N)$. The random Fourier features (RFF) technique, based on sampling from the spectral density of the kernel function, has become popular to speed up kernel estimators for machine learning applications. Unfortunately, it is restricted to positive definite kernels, while the majority of kernel functions popular in KDE, such as the parabolic kernel, do not satisfy this property. To overcome this limitation, this article introduces the signed random Fourier features (SRFF) technique. It is a generalization of RFF compatible with indefinite kernels whose inverse Fourier transform is absolutely integrable. The motivation for introducing this method is to speed up KDE in the case of multivariate compact kernels, which are generally not positive definite. We detail how to implement SRFF for both product kernels and isotropic kernels. For the class of Kuttner-Golubov kernels $K(\boldsymbol{x}_{i},\boldsymbol{x}_{j})=(1-\left\Vert \boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\Vert ^{\alpha})^{\beta}\mathbf{1}_{\{\left\Vert \boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\Vert \leq1\}}$ where $\boldsymbol{x}_{i}\in\mathbb{R}^{d}$, $\boldsymbol{x}_{j}\in\mathbb{R}^{d}$, $\alpha>0$, $\beta>0$, which includes the triangular, parabolic, biweight, triweight, and other kernel functions of interest for KDE as particular examples, we provide an explicit acceptance-rejection algorithm to sample from its signed spectral density. Our numerical tests on a dataset of one million points confirm the computational efficiency and accuracy of SRFF for large-scale KDE.

Varmapack is a software package for efficient, exact simulation of VARMA time series without a burn-in period. For stationary models, Varmapack can generate initial states and innovations from their joint stationary distribution. Alternatively, the user can supply initial states, in which case innovations are generated from their conditional distribution. The core C library has interfaces to R, Python, and Matlab. It also offers VARMAX simulation and computation of autocovariances and correlations, spectral radii, and impulse responses. Varmapack automatically selects between vector-Yule-Walker and state-space methods for covariance computation and uses level-3 BLAS operations for efficient generation of multiple replicates. Benchmarks on several platforms show substantial performance gains over existing simulation software, ranging from several-fold to more than three orders of magnitude for the packages and models considered. Varmapack is open source and publicly available through GitHub.

We study the mixing time of weighted Dikin walks for sampling from exponential distributions on polytopes and truncated positive-semidefinite (PSD) cones. Our first result gives a general total-variation mixing bound under strong self-concordance, $\bar{\nu}$-symmetry, and mixed-trace regularity on the local metric. The key idea is to control the Metropolis–Hastings acceptance probability on a high-probability region rather than at every point. Applying this framework to the Lee–Sidford, Lewis-weight, and John metrics yields an $\widetilde O(d^2)$ mixing bound for sampling from polytopes, while applying it to a hybrid barrier yields an $\widetilde O(d^4)$ mixing bound for sampling from truncated PSD cones. Our second result establishes stronger $\chi^2$-divergence guarantees and pointwise acceptance control using a new fourth-order bootstrap condition. For a suitably scaled Lee–Sidford metric, this yields an $\widetilde O(d^2)$ mixing bound in $\chi^2$-divergence, improving on the previous $\widetilde O(d^{9/4})$ bound.

Relational Event Models (REMs) provide valuable insights into the dynamics of longitudinal social networks. Yet, the vast availability of potential predictors for a dyad's event rate poses the risk of selecting irrelevant variables and specifying an overfitted model that does not generalize to new data. Despite the recent popularity of Bayesian regularization methods to address this, there has not been a systematic evaluation of the performance of Exact Bayesian Regularization (EBR) and Approximate Bayesian Regularization (ABR) against standard Maximum Likelihood (ML) procedures commonly used for REMs. To address this, we conduct a simulation study in which directed Relational Event History (REH) data is generated with endogenous and exogenous effects of varying strength and compare the performance of (a) unregularized ML estimation, (b) ABR through normal approximations of the likelihood with Ridge and Horseshoe priors, and (c) EBR with Horseshoe priors. We find that threshold-based criteria applied to ABR and EBR with Horseshoe priors outperform univariate selection using p-values from ML estimation in variable selection accuracy. Regarding predictive performance, ABR with Ridge priors and EBR with Horseshoe priors outperform ML REMs, particularly for small sample sizes. Given only small gains with large computational burdens when using EBR, we therefore advise researchers to use ABR with a suitable prior.

Modern time series are often long, serially dependent, and non-stationary. Existing change-point methods either target specific changes or become computationally intensive when using nonparametric costs on long series. Many also require thresholds to be carefully calibrated under serial dependence. We introduce SCAN, an offline method for detecting multiple distributional change-points in long, serially dependent univariate time series. SCAN compares adjacent windows using an integral probability metric, calibrates local discrepancies with a dependence-aware bootstrap, and refines candidate locations using a scaled 1-Wasserstein criterion, enabling detection of changes in mean, variance, and broader distributional structure within a unified framework. An ensemble over multiple window sizes reduces sensitivity to window size and threshold specification. We establish consistency of the estimated number and locations of change-points under exponential alpha-mixing dependence, and show that the localization statistic reduces to a CUSUM-type statistic under pure mean shifts. In simulations with up to one million observations, SCAN generally achieves higher covering and F1-scores than competing methods across mean and joint mean-variance shifts, particularly under serial dependence. On real data, SCAN identifies labeled activity transitions in sensor data and interpretable structural changes in hourly Bitcoin prices. Implementations are available in the Python package scan-cpd and R package scanr.

We present a suite of R packages for macroeconomic forecasting that leverages advanced Bayesian, structural, multivariate, dynamic, hierarchical, non-linear, and non-Gaussian models. The suite enables both structural and predictive analyses, and is adapted to time series data across various types, dimensions, and sampling frequencies. Each additional feature increases computational complexity. To address this challenge, our software design incorporates a carefully curated selection of models, efficient algorithms implemented in C++, advanced econometric and numerical methods, robust handling of complex input and output objects, and standardised workflows. This approach combines the computational efficiency of C++ with the convenience of working with data in R. We demonstrate that our packages facilitate original research contributions in forecasting, as illustrated by our example in which vector autoregressions with non-centred stochastic volatility enhance density and point predictions relative to models with centred stochastic volatility.

The problem of sampling from a probability distribution arises in many applications such as posterior sampling in hierarchical Bayesian inverse problems and Gaussian processes for machine learning. Markov chain Monte Carlo (MCMC) algorithms are often used for sampling from a target probability distribution, but implementations can be computationally expensive, especially for large-scale problems. In certain applications, the target distribution naturally factorizes into a lower dimensional marginal distribution and a conditional distribution that allows exact sampling. We describe an MCMC algorithm called MarCo that exploits such a structure and generates a Markov chain via Metropolis-Hastings sampling from the marginal distribution, followed by sampling from the exact conditional distribution. By design, MarCo constructs a Markov chain on the joint space that inherits the convergence behavior of the marginal MCMC algorithm. This provides multiple theoretical and computational advantages. We prove that MarCo can achieve improved mixing time upper bounds compared to direct sampling from the joint distribution. Moreover, compared to one-block methods that also exploit marginal-conditional structure, we use the framework of Peskun-Tierney ordering to show that MarCo has a larger right spectral gap and smaller asymptotic variance, thus leading to superior convergence properties. Numerical results illustrate the performance benefits of MarCo and are provided for various problems, including a semi-blind image deblurring example.

High-dimensional semilinear parabolic partial differential equations arise in stochastic control, financial engineering, and uncertainty quantification, but classical spatial discretizations suffer from the curse of dimensionality. Motivated by Gaussian perturbation and conditional regression in denoising score matching, we propose Deep-Control BSDE (\DCBSDE{}), a layerwise control-regression method for Markovian backward stochastic differential equations. From an implicit Euler scheme, we identify the discrete ideal control $z_n^\pi=\mathcal M_nu_{n+1}^\pi$ as the conditional projection coefficient of the successor value response onto the one-step Brownian increment. At each time level, the method freezes the successor value function, approximates the resulting Brownian conditional-moment target using finitely many branches, regresses the control, and then fits the value through the implicit BSDE relation. A frozen linear-response baseline and antithetic pairing reduce finite-branch fluctuations while preserving the conditional target, and a within-layer correction coordinates the value and control approximations. We establish backward stability through a contraction property of the Brownian projection and derive conditional consistency when discretization, local learning, numerical, and finite-branch errors vanish jointly. Experiments across six benchmarks demonstrate that \DCBSDE{} achieves a favorable overall balance among value accuracy, control accuracy, and dynamic consistency while exhibiting generally stable performance across random seeds. These results support the use of \DCBSDE{} in applications requiring reliable joint approximation of the value and control processes.

Joint analyses across multiple institutions are increasingly important in biomedical and epidemiological research, particularly for rare diseases where datasets are typical small. However, privacy regulations and institutional policies often prevent the sharing of individual-level patient data. In this paper we present an accurate and single-communication federated inference algorithm. Single-communication federated inference enables statistical analyses through a single exchange of summary statistics between participating centers and a coordinating server, preserving privacy while reducing communication and computational costs compared with iterative federated learning. We extend a recently proposed single-communication federated inference strategy that is based on second-order Taylor expansions by using third-order expansions to better approximate local log-likelihood functions. The proposed method is evaluated through simulation studies based on real data and compared with existing federated inference strategies. The simulation studies assess the performance of the proposed method, with a particular focus on scenarios involving small local sample sizes, where quadratic approximations may fail to capture skewness and other higher-order characteristics of the log-likelihood function. They demonstrate that incorporating higher-order information of the log-likelihood function improves the accuracy while preserving the privacy, communication efficiency, and scalability required for collaborative biomedical and epidemiological research.

Goodness-of-fit tests for logistic regression are routine in clinical risk modelling, yet the classical tests say whether a model misfits, not where. A pretrained network is not a test until its level, power and blind spots are established. We audit DeepGOF-1, which renders the residuals of a fitted logistic model as a map over the ranks of two covariates, scores the map with a frozen convolutional network, and calibrates the score by the analyst's parametric bootstrap. The level is first-order valid under standard bootstrap regularity and two conditions, consistency against a named alternative is decided by one forward pass, and local power has an explicit limit. In simulation the level holds from 20 patients to 8,873; the network adds power over a chi-square statistic of the same map when misfit is sparse among many covariates; the map finds a missed interaction that one-covariate diagnostics cannot see; and, applied to a published risk model on new patients, it is an exactly valid test of calibration within patient subgroups, where the calibration belt has little power. In the SUPPORT study of 8,873 inpatients, DeepGOF-1 rejects a linear in-hospital mortality model, its map shows the U-shaped risks of blood pressure and respiratory rate that a calibration curve hides, and it guides the repairs. It ran in 13 seconds; the projection test, the strongest omnibus rival, took 9 hours 53 minutes and 4.5 GB on the same data, too slow for routine use. The test is deepgof1() in the R package ebrahim.gof.

Stein's method constructs computable discrepancies between a target distribution and a candidate distribution without requiring the target distribution's normalizing constant. These discrepancies support goodness-of-fit tests for model assessment as well as sampling tools for empirical approximation. The R package steinsampling provides the first unified R workflow for applying score-based Stein methods to kernel goodness-of-fit testing of independent or serially dependent observations, point transport, greedy point construction, and sample compression. High-level functions carry out each task in a single call, while the kernel, calibration, optimization, and transition components are provided separately so that users can replace any one of them. A single score and kernel setup can therefore be reused across sampling and testing, making these methods easier to reproduce, compare, and extend.

The OBABO and BAOAB schemes and the other standard Strang splittings of kinetic (underdamped) Langevin dynamics are widely used Markov chain Monte Carlo algorithms. Under a suitable friction scaling, the underlying diffusion relaxes on a ballistic time scale, suggesting that these discretizations, suitably tuned, sample targets with condition number $\kappa$ in $O(\sqrt{\kappa})$ iterations. We prove that no fixed choice of step size and friction, based only on the curvature bounds and the dimension, achieves this acceleration: total variation mixing time lower bounds for OBABO show that ballistic cold-start mixing fails uniformly over the smooth strongly convex class, and the lower bounds extend, with the same orders, to BAOAB and the other four Strang splittings. The proof transfers non-acceleration from optimization to sampling. Eliminating velocity gives an exact noisy heavy-ball recursion, and by the non-acceleration theorem of Goujaud, Taylor and Dieuleveut, for every tuning either some Gaussian target has a mode with relaxation time at least of order $\kappa$, or an attracting cycle exists on a smooth potential; dilating such a potential as $U_R(x)=R^2U(x/R)$ preserves its curvature bounds and produces metastability for a number of steps exponential in $R^2$, from an initial state at Wasserstein distance $O(R)$ from equilibrium. Using contraction estimates of Leimkuhler, Paulin and Whalley and a Wasserstein-to-total-variation regularization estimate, we prove a complementary upper bound of $O(\kappa)$ steps, up to logarithmic factors, for a fixed-parameter OBABO tuning. Hence, among fixed-parameter OBABO tunings, the optimal condition-number dependence of cold-start total variation mixing over this class is linear, up to logarithmic factors. A direct Gaussian calculation also rules out fixed-parameter acceleration for the left-endpoint exponential integrator.