Papers

We introduce a sampling approach for energy- and score-based generative models that requires no gradient evaluations of the model. Replacing the drift term that would normally contain the score $\nabla_\mathbf{x} \log p_\theta(\bf{x})$ with a high-frequency dithered cosine of the model's value, $\sqrt{\alpha\omega}\,\cos(\omega t + k \log p_\theta(\bf{x}))$, produces, in the high-frequency averaging limit, Langevin Markov chain Monte Carlo for energy-based models and the reverse-time SDE of score-based diffusion. We prove that trajectories of the dithered Itô SDE converge to those of the target SDE, driven by the same Brownian motion, uniformly on compact time intervals in probability, by an averaging argument that extends bounded extremum seeking (ES) to Itô processes, with an explicit $O(\omega^{-1/2})$ mean-square rate under global bounds. The approach is not confined to smooth targets: it extends to $C^{1,1}$ energies with discontinuous curvature (without ellipticity requirement) and to Sobolev energies whose Hessians exist only off measure zero sets; for Lipschitz energies with gradient kinks the averaged limit remains well posed; the Krylov-Röckner integrability class is the boundary of provability. The approach provides a hard a priori bound on the per-step update rate and applies to explicitly time-varying targets on finite horizons. Gradient-free pixel-space sampling is not competitive with well-tuned backpropagation-based samplers at practical evaluation budgets; the regime where the approach offers an advantage is latent-space sampling when the model is a black box and the target drifts in time. We demonstrate latent-space tracking for time-varying images on CelebA-HQ ($256{\times}256$) from limited 1D projection measurements and latent-space EBM sampling on CIFAR-10.

Fully homomorphic encryption allows arithmetic to be carried out on encrypted values, so that a party performing a computation need not see the data it operates on. This is useful wherever collaborating parties or sites cannot share data records or computed summaries yet need to do a joint analysis. We demonstrate the protocols needed to perform such analyses reproducibly via two packages in R, a platform widely used for applied statistics. The first, openfhe.R, is an interface to the OpenFHE C++ library, which exposes exact integer arithmetic (BFV, BGV), approximate real-valued arithmetic (CKKS), Boolean circuits, and multiparty key generation. The second is homomorpheR, which builds a small set of master/worker primitives for multi-party protocols on top of the first. Together, the two let an ordinary R statistical routine compose with an encrypted-arithmetic layer at the function-value boundary, provided the quantity crossing that boundary decomposes as a sum over sites. Two validation studies are presented: distributed stratified Cox regression, where stats4::mle() converges through an encrypted channel, and federated Cox-lasso via convex optimization using consensus ADMM, where the cross-site update is the only encrypted step. We also discuss federated similarity retrieval and two-party prediction. Throughout, the parties are assumed honest-but-curious, and threshold key generation removes the assumption that any single party holds a usable secret key.

The LSJM R package provides a series of functions to estimate models based on location-scale linear mixed model and location-scale joint model theory. For longitudinal data, it includes the estimation of location-scale mixed models with subject-specific residual variance that can be time-dependent or distinguish within-visit and between-visits variability. When survival data is considered, implemented joint models handle competing risks as well as semi-competing risks possibly with interval-censored data through the use of a illness-death sub-model. They incorporate flexible dependence structure between the longitudinal marker and the times-to-events through the current value, the slope and/or the residual variability of the marker. Maximum likelihood estimators are obtained using the Marquardt-Levenberg algorithm with stringent convergence criteria. The package also provides several post-fit functions including goodness-of-fit analyses and individual dynamic predictions. Some examples on real data allows to illustrates the use of this package.

With the advance of measurement technologies and increasing computing power, large spatial data with heterogeneous structures are often collected over high-dimensional domains. Existing Gaussian process (GP) models and computational strategies are often inadequate for analyzing such datasets in multi-dimensional domains. To address these challenges, we develop a Bayesian residual tree GP methodology called ResTGP for large spatial data with potentially heterogeneous structures in multi-dimensional domains. The key idea is to decompose a Gaussian process at a cascade of resolutions along a dyadic tree through iteratively computing predictive and residual processes so that the residual process on each tree node, both interior and leaf, becomes sufficient for the finer-level dependency within that node. This allows characterization of the underlying covariance structure in a flexible, multi-scale manner while achieving divide-and-conquer on the data domain, which leads to computational efficiency. To allow efficient tree inference, we introduce a computational strategy for Bayesian inference based on recursive message passing, which scales linearly with the sample size given the tree. This paper also proves posterior consistency of the model for estimating continuous functions in a nonparametric regression framework. Extensive numerical examples and the storm surge application confirm the advantages of the proposed method.

Sampling from the posterior is the central task of computational Bayesian inverse problems. The standard workhorse in Bayesian inference - Markov chain Monte Carlo (MCMC) - is sequential, yields correlated samples, and must be rerun for each observation. Conditional flow matching offers an amortized alternative: a transport map, trained once on joint samples of parameter and data, that yields independent approximate posterior samples for any observation at negligible online cost, without new likelihood evaluations. We give a careful, MCMC-literate assessment of flow matching for PDE-based inverse problems with function-valued parameters. Exploiting the flow's tractable density, we derive computable accuracy estimates of the underlying approximate posterior in total-variation distance and Kullback-Leibler divergence and, moreover, propose a hybrid sampler that is asymptotically exact by Metropolization. We validate the accuracy estimates and demonstrate the amortization in several numerical examples, including electrical impedance tomography and a likelihood-free Lotka-Volterra model.

We introduce a Monte Carlo approach for estimating isolated polynomial zeros through their contour integral representations. For a simple zero enclosed by an isolating contour, we show that the zero can be expressed exactly as the expectation of a complex-valued random variable obtained by uniformly sampling the contour. This yields an unbiased estimator without perturbing the polynomial coefficients. We establish its variance and derive finite sample probabilistic error bounds based on concentration inequalities, with the bounds explicitly reflecting the geometry of the contour and its separation from the zeros. The framework also provides a stochastic formulation of contour based root counting and root isolation. Numerical results illustrate the estimator's behavior under varying polynomial degrees, and we further provide comparisons to existing methods.

Differential privacy implementations rely on precise sampling from noise distributions to provide formal privacy guarantees. We report the discovery of systematic artifacts in OpenDP's discrete Laplace sampler that manifest as periodic distortions in the output distribution. Through systematic testing, we trace these artifacts to a faulty implementation in the rational arithmetic library used by the bernoulli_exp1 function, a low-level primitive that implements sampling from Bernoulli(e^(-x)) distributions. We present a diagnostic methodology that isolates the faulty component in the nested sampling hierarchy and propose an alternative implementation based on exact rational arithmetic that eliminates the artifacts. Statistical validation with 10^6 samples confirms that the corrected sampler produces outputs indistinguishable from the theoretical distribution at the tested precision level.

Influential Outlier Detection is developed for mixed-effects models on clustered data. The Influential Outlier Metric is defined as a combination of SHapley Additive exPlanantion (SHAP) values and model residuals, both of which undergo a change of measure transformation. Building on previous work showcasing the suitability of using Normalizing flows to map arbitrary distributions to a flexible base distribution for statistical inference, the Normalizing Flows are constructed to allows contextual information and also provide a goodness of fit diagnostic for model evaluation. The use of SHAP values in the construction moves away from model specific tools and instead provides point-wise model agnostic influential outlier. The advantages and limitations of this approach are examined in several models including the linear model, the random forest, and gradient-boosted trees.

Event-chain Monte Carlo (ECMC) has revolutionised computational sampling over recent years, providing a powerful alternative to the molecular-dynamics (MD) and Hamiltonian Monte Carlo (HMC) algorithms. Each method outperforms the ubiquitous random-walk Metropolis algorithm by advancing particles along deterministic trajectories, but ECMC achieves this without being constrained by Newtonian dynamics. Recent advances exploited this dynamical freedom to induce a collective particle dynamics that relax local density variations on fast timescales. In a foundational model of $N$ pairwise-interacting particles on the 1D torus, this led to an $O(N^{3/4})$ improvement on the industry-leading computational efficiency of MD and HMC - but with an impractically small prefactor at low/high mean particle density for general attractive/repulsive interactions. Here we present a universal framework that generalises this high-efficiency sampling strategy to all translationally symmetric pairwise models on the torus - creating the potential to surpass MD and HMC as the state of the art. We also numerically elucidate the collective particle dynamics and discuss broad impact across the physical sciences and computational statistics.

We prove asymptotic bias bounds for stochastic gradient Langevin dynamics in Wasserstein distance of order two. We assume that the negative log-density is strongly convex with a Lipschitz gradient, and that the stochastic gradient estimator is unbiased with an error satisfying a mean-square Lipschitz condition. The bounds are of order $h$ under a fourth moment assumption on the stochastic gradient error and of order $h^{1/2}$ under only a second moment assumption, where $h$ is the stepsize. A spiked-noise example shows that a second moment assumption alone is insufficient for a bound of order $h$ that is uniform over noise distributions with a fixed variance.

Modern machine-learning models often contain strongly dependent or redundant features, making feature attribution difficult because shared predictive information can be distributed across correlated predictors. Existing methods such as SHAP, LIME, HSIC, MI/CMI, and SAGE may therefore produce unstable rankings under multicollinearity or near-duplicate predictors. We propose the Mutual Correlation Impact Ratio Method (MCIR-M), a dependence-aware global feature-importance approach that quantifies the unique predictive information contributed by each feature beyond a selected dependence neighbourhood. MCIR-M introduces the Mutual Correlation Impact Ratio (MCIR), which conditions each feature on strongly dependent neighbours and computes a normalized ratio of conditional to block-level information. The population score lies in [0,1] and equals zero under exact conditional redundancy. We also introduce a lightweight estimation procedure that computes MCIR using a fraction of the available data and evaluates agreement with full-data explanations. Across controlled synthetic redundancy experiments and the UCI HAR benchmark, MCIR shows dependence-aware ranking behaviour, with its clearest advantage under injected near-duplicate predictors. Comparisons with independent and conditional SHAP, SAGE, HSIC, MI-based scores, and CIR-family baselines are mixed across real-data criteria. Reduced explanation samples lower computational burden in the evaluated configurations, while agreement with full-data explanations is assessed separately through ranking, head-set, and faithfulness diagnostics. Overall, MCIR-M provides a practical dependence-aware diagnostic for global explanation under strong feature dependence.

Skew-symmetric probability distributions provide a principled mechanism for incorporating gradient information into Markov chain Monte Carlo algorithms. Here we review the (preconditioned) Barker proposal, a Metropolis–Hastings algorithm built on skew-symmetric distributions, and motivate its design. We then introduce three natural extensions. First, we propose coordinate-free variants of the Barker algorithm. Second, we introduce a Gibbs-style Barker algorithm that re-evaluates the gradient at each partially updated coordinate. Third, we derive a simplified manifold Barker algorithm, producing a manifold sampler with enhanced robustness compared to natural comparators. Numerical experiments demonstrate that the Gibbs-style variant improves raw sampling efficiency on correlated targets, that the coordinate-free variants offer limited practical advantage over the standard Barker proposal once computational costs are accounted for, and that the simplified manifold Barker algorithm can achieve significant advantages over simplified manifold MALA when the local geometric structure of the target is irregular or unreliable.

Gibbs samplers, and more generally component-wise Markov chain Monte Carlo algorithms such as Metropolis-within-Gibbs, can be implemented using either random-scan or deterministic-scan updates. How much convergence can depend on this choice of scanning rule has been a longstanding question. We study this problem through $L^2$ spectral gaps, measuring computational cost in units of component updates. For a $d$-block Gibbs sampler, the cost of random scan is at most twice that of any fixed deterministic scan, while the reverse cost ratio is at most of order $d^2$. We extend these comparisons to general reversible component-wise updates under the global block-wise contraction condition. If $K_j$ denotes the update of block $j$ and $P_j$ its Gibbs counterpart, and $\|K_j-P_j\|\leq \lambda_0<1$ for $j=1,\ldots,d$, then the cost of random scan is at most $2/(1-\lambda_0)$ times that of deterministic scan, while the reverse cost ratio is of order at most $d^2/(1-\lambda_0)$. Examples show that the cost bounds for random scan relative to deterministic scan are asymptotically sharp and that the joint dependence on $d$ and $(1-\lambda_0)^{-1}$ in the reverse comparison cannot be improved uniformly. This work was assisted by generative AI, including for formal verification of mathematical results in Lean. The human author reviewed and verified the mathematical content and takes full responsibility for the results.

Speculative sampling accelerates diffusion generation by verifying inexpensive draft states in parallel while preserving the target law. Recent tree-based methods allocate the parallel compute budget more effectively than single-chain drafts, as demonstrated by Diffusion Greedy Rejection Sampling (D-GRS). D-GRS generates $K$ conditionally independent candidates per node, and sequentially tests them in their generation order. Yet the sampled candidates admit an informative ranking without additional target-model evaluations. To exploit this, we introduce Rank-Aware Speculative Sampling (RASS), a verification rule for speculative draft trees based on rank-aware list coupling. RASS orders draft candidates along the proposal-target mean displacement and samples a rank with weights optimized to minimize total variation between the selected-proposal and target laws. Finally, the selected candidate is maximally coupled with the target, with residual correction ensuring exact sampling for any choice of rank weights. We evaluate RASS on a Gaussian-mixture target, unconditional pixel-space generation on FFHQ, conditional generation on CIFAR-10, and latent diffusion with Stable Diffusion 3.5 using COCO2014 prompts. Measured by the ratio of standard to speculative sampling's target-model evaluation counts, RASS improves on D-GRS across the evaluated settings, with gains reaching approximately 20% on CIFAR-10 at matched compute budgets.

We develop multifidelity methods for constructing triangular transport maps from samples, when high-fidelity data are scarce but lower-fidelity data are more abundant. Using this set of multifidelity data, we approximate a triangular transport map that bijectively maps between a tractable reference density and the high-fidelity target distribution. We introduce two strategies to leverage low-fidelity data: a hierarchical approach that composes maps between adjacent fidelity levels, and a non-hierarchical method that incorporates low-fidelity information through monotonicity-preserving corrections to the map parameterization. Numerical experiments compare these strategies with single-fidelity transport and demonstrate how the proposed multifidelity approaches can improve map estimation from limited high-fidelity data. To illustrate the broader utility of the learned maps, we also deploy them in a downstream amortized simulation-based inference task. This example shows that multifidelity improvements in map estimation can translate to improved conditional sampling and uncertainty quantification when high-fidelity data are scarce.

We develop a general method for amortized Bayesian inference on multilevel models of arbitrary structure. Given a generative model specified as a directed acyclic graph, our method automatically derives valid factorizations of the joint posterior and matching neural network architectures. The key steps, graph expansion and graph inversion, yield an inverse graph that determines how inference networks are stacked and conditioned, producing factorizations that amortize over the number of groups and the number of observations within each group. Unlike approaches that simplify the dependency structure to speed up learning or inference, our method preserves all conditional independence and exchangeability assumptions of the generative model. Across three case studies, it closely matches gold-standard samplers on models with more than 6,500 parameters while reducing inference to a near-instant forward pass once trained.

We consider simulation-based Bayesian inference (SBI) for the parameters of models with intractable likelihoods but tractable forward simulation. Building on Gaussian mixtures-of-experts surrogates, as a first contribution we develop a sequential procedure to estimate the posterior, incorporating an importance sampling step that introduces a locally amortized likelihood-to-evidence ratio estimator. The local amortization allows the ratio estimator to be evaluated at new parameter–data pairs from a local (data-dependent) design, without re-estimation, which is particularly convenient for our second contribution. The second contribution is localized simulation-based calibration (SBC). Localized SBC provides a graphical diagnostic of calibration over neighborhoods that are deliberately broader than the available posterior, at low additional computational cost. Rather than repeatedly running the full inference procedure for each simulated pseudo-observation, localized SBC fits the local surrogate and ratio estimator only once for each local design, and reuses them across all pseudo-observations. This approach is particularly appealing for SBI applications. We evaluate both contributions on three case studies, including a real-data epidemiological application.

Casting Bayesian inference as a neural network optimization problem targeting an amortized posterior is attractive, as it extends to otherwise intractable statistical models and offers near instantaneous inference for new datasets after prepaying the training cost. Although theory guarantees faithfulness under ideal convergence, practical amortized inference still requires iterating over architectures and optimization choices and ultimately “satisficing” under finite simulation, compute, and time budgets. Even the best-performing solution may thus retain avoidable representation gaps that typically require problem-specific fixes. Here, we propose a generic alternative which improves training dynamics with auxiliary guidance losses applied to internal representations. Specifically, we show how such guidance leads to faster convergence when training data is abundant and to better performance when it is scarce. We formalize representation gaps as getting stuck in a local optimum at the information bottleneck between the parts of the network tasked with feature learning and those tasked with conditional distribution learning, and offer a generic diagnostic to separate summary failures from inference failures. Finally, we demonstrate that auxiliary supervision improves convergence speed and accuracy on a range of challenging real-world inference problems.

We present B-GRASP (Bayesian GRAph inference with SPDE priors), a Bayesian framework for inferring uncertain edge weights in stochastic dynamical systems on graphs from noisy observations of nodal states. The unknown edge weights parameterize the graph differential operator and therefore directly govern the evolution of the graph process. Motivated by connections between differential operators in PDEs and SPDEs and their graph counterparts, we construct stochastic graph models incorporating diffusion, reaction dynamics, and stochastic forcing. The resulting hierarchical formulation jointly represents uncertainty in the graph structure and stochastic forcing. Latent graph variables determine positive edge weights and the corresponding graph Laplacian, while latent Brownian variables represent the stochastic forcing. Conditional on these variables, the graph dynamics define a deterministic forward map from which the likelihood and posterior distribution are constructed. We characterize the posterior using maximum a posteriori estimation and the No-U-Turn Sampler, enabling both point estimation and uncertainty quantification. We demonstrate the framework on a one-dimensional inverse heat-conduction problem, stationary and nonstationary graph reaction-diffusion systems with nonlinear dynamics, and state-level COVID-19 data in the United States. The numerical results show that posterior uncertainty provides information not captured by point estimates, particularly for weakly identifiable or highly conductive edges, and enables uncertainty in both graph connectivity and stochastic forcing to be quantified within a Bayesian framework.

Learn-then-differentiate (LTD) estimates gradients by fitting a model to simulation outputs and differentiating it. We develop a unified framework explaining what LTD differentiates and how accurately it estimates gradients. For models with a weighted representation, LTD differentiates a learned representation of the underlying probability measure. We then show how accuracy guarantees for fitted models translate into guarantees for gradients and higher-order derivatives, with rates approaching the standard Monte Carlo rate under suitable smoothness conditions. The framework recovers established results for kernel regression, local polynomial regression, and kernel ridge regression, and yields further guarantees for multiple kernel learning and smooth neural networks. These results provide a common foundation for understanding and analyzing LTD across learning methods.

Structured additive distributional regression flexibly relates all parameters of a conditional response distribution to covariates, but multivariate extensions remain challenging when dependence is complex. We propose a multivariate structured additive distributional regression model based on regular vine copulas. Different vine edges may use different pair-copula families, while every pair-copula parameter may vary with covariates through a structured additive predictor. The model therefore accommodates heterogeneous marginal distributions together with pair-specific asymmetric, tail-dependent, and covariate-dependent dependence structures. For scalable inference, we develop a tree-wise stochastic variational inference procedure based on component-specific Gaussian variational approximations. Marginal models are estimated first, followed by pair-copula regressions sequentially along the vine trees. We also adapt sequential vine selection by fitting candidate covariate-dependent pair-copula regressions and using information criteria based on effective degrees of freedom for both family and tree selection. In simulations, the tree-wise estimator remains close to an oracle using the true recursive conditional inputs, whereas global refinement yields more concentrated approximations and lower frequentist coverage for upstream components. An application to six-dimensional meteorological data from the Netherlands yields a vine combining Gaussian and non-Gaussian pair copulas, with pronounced nonlinear spatial and temporal variation in dependence.

Multimodality is common in many scientific and engineering applications. However, standard Markov chain Monte Carlo (MCMC) convergence diagnostics often penalise multimodality in practice by conflating overall poor mixing with groups of chains that have mixed well within distinct modes. To diagnose multimodality rather than flag it as a sampling failure, we introduce the MCMC convergence graph: a graph-based diagnostic that computes pairwise R-hat values between chains and summarises them in a graph, revealing sets of chains that explore the same posterior mode. We implement the method in the probabilistic programming framework Stan and demonstrate its use on worked examples ranging from regression models with inherent multimodality in the data, to a pharmacokinetic model for which multimodality reflects non-identifiability. The graph-based diagnostic, which we make freely available in the R package mcmcConvergenceGraph and which accepts MCMC output from any sampler, provides both graphical and numerical summaries of MCMC behaviour for diagnosing multimodality.

The paper introduces the R package CovEsts, which implements several nonparametric estimators for the autocovariance function. First, it presents the theoretical foundations of the implemented estimators, their properties, assumptions, and potential limitations. Next, it outlines the structure of the package and its key functions, including several methods for estimating autocovariance functions, constructing corresponding bootstrap confidence regions, and correcting the provided estimators. The package also includes diagnostic tools, such as several metrics for comparing estimators, and additional functions for broader use. The article illustrates a high degree of flexibility of the package in the selection of function parameters and the tuning of the estimators. Applications of selected estimators and package functions are illustrated using simulated data, yearly sunspot counts, and US unemployment increments data.

In Bayesian inference problems with non-Gaussian observation noise, the posterior is only as accurate as the noise density, and gradient-based samplers need that density and its gradient evaluable pointwise, whether from an explicit expression or from code, and without an inner solve. We propose Copula Active Subspaces (CAS) to represent this noise density. A componentwise rank transform isolates the noise law's dependence in its copula, and a rank-$r$ reduction keeps only the directions along which that dependence varies. These directions are the leading eigenvectors of the copula score covariance $\boldsymbol{C} := \mathrm{Cov}_{\pi_{\boldsymbol{Z}}}(\nabla\log c^{Z})$, which is what makes the reduction a copula active subspace. Because $\boldsymbol{C}$ vanishes when the coordinates are independent, these are directions of dependence, which the covariance of the data need not identify. From this construction follow a Gaussian-reference KL divergence bound with the explicit constant $\tfrac{1}{2}$, minimized over all rank-$r$ reductions by exactly this eigenspace; a diagnostic for the error the reduction leaves behind, computable from the samples alone; and, from Hermite score matching, a reduced log-density and gradient in closed form, with the truncation orders and the Stage-2 regularization constants chosen on validation samples. The reduction replaces a $d$-dimensional density estimation problem by an $r$-dimensional one. On a $d=20$ noise law and a Bayesian inference problem with that noise, CAS lowers noise KL divergence more than fivefold and posterior KL divergence more than sevenfold against Gaussian-copula, product-of-marginals, and PCA-subspace baselines, and lowers noise KL divergence by factors of about $3.5$ and $2.7$ on two further $d=20$ examples.

Many inference-time tasks for pretrained discrete diffusion models and diffusion language models reduce to drawing samples from a tilted version of the pretrained distribution. Feynman-Kac sequential Monte Carlo (SMC) makes this correction exact in principle, but its prescribed weights routinely degenerate when the proposal dynamics are misaligned with the tilt, capping the practical gains from additional particles. We introduce FluxLite, a lightweight, training-free proposal-control framework for discrete diffusion. On the sparse directed graph of pretrained reverse rates, any sparse jump-rate perturbation can be exactly compensated by a $q_t$-weighted graph-divergence term in the Feynman-Kac potential; the target path is therefore preserved while the residual reweighting variance becomes a local convex objective. We instantiate this principle as two practical samplers: a one-hop local reallocation rule (HEU) and a small nonnegative quadratic program over pretrained-rate bases (D-VCG). We further prove population stability under the standard score-entropy training loss, identifying a tilted-path coverage factor that governs robustness to score error, together with finite-particle convergence for a fixed controlled Feynman-Kac recursion. Empirically, FluxLite improves over standard Feynman-Kac SMC baselines by up to two orders of magnitude in terminal KL on an analytically tractable finite-state CTMC benchmark, and reduces row-correlation MSE on 2D Ising sampling by 5-7x in geometric mean and up to 55x at peak.

Response-adaptive randomization (RAR) modifies treatment allocation probabilities during a clinical trial as response/outcome data accumulate, with the aim of improving patient benefit, statistical efficiency, or both. Despite substantial methodological development, adoption of RAR in clinical practice has remained limited, and the software available to support its design and implementation has not previously been reviewed. We identified 16 publicly available, open-source software packages implementing RAR methods, spanning urn-based, target-allocation, Bayesian, Markov decision process (MDP)-based, and dose-finding approaches, and evaluated them with respect to their methodological and practical characteristics. We found that while several software packages exist, most are method-specific, and only a small number provide broader, general-purpose adaptive-trial design and analysis capabilities. Explicit support for practically relevant features, including delayed and missing outcome data, temporal trends in response rates, flexible operating-characteristic evaluation, and platform or multi-arm multi-stage trial designs, is rare or absent across the identified software. These findings suggest that while a diverse set of tools exists for exploring RAR designs, gaps remain between the methodological literature and the software available to implement it in practice. Continued development of flexible, practically oriented, and validated software is important for the wider adoption of RAR in clinical research, and we highlight interesting areas of further work.

In ecology, psychometrics, and the analysis of social and financial networks, binary matrices are often analyzed conditional on their observed row and column sums, which restricts the problem to a finite sample space of matrices with the same margins. Two fundamental problems are to count this space and to sample uniformly from it. Sequential importance sampling (SIS) addresses both with independent weighted samples and an unbiased count estimator, but its efficiency depends critically on the proposal distribution. Existing proposals are analytically designed, and their accuracy can vary substantially with the margins. We show that the ideal SIS proposal, under which every weight equals the count and the variance vanishes, is exactly the policy of a generative flow network (GFlowNet) with unit reward on every matrix that has the given margins. We therefore propose MarginFlow, a framework that turns the design of the proposal into a learning problem and amortizes it across margins by exploiting their self-similarity. Every partial matrix is itself an instance with reduced margins, so one set transformer that reads the remaining margins serves every margin. We train MarginFlow on a pool of 1904 margins and evaluate it zero-shot on 1190 held-out margins, synthetic and real, from $3\times3$ to $870\times6$. On 1187 of the 1190 margins it matches or beats the best of 31 analytically designed configurations, chosen post hoc for each margin, and its median effective sample fraction is 99.8%. On the 56 margins where that best loses more than one nat of effective sample size, MarginFlow wins every one and raises the median effective sample fraction from 10.3% to 94.1%.

Virtual simulators are widely used for studying complex physical phenomena, from particle collisions to rocket propulsion. Such "computer experiments" can be highly time-intensive, and a Bayesian surrogate model can be used for efficient emulation with reliable uncertainty quantification. To train accurate surrogates with a limited sample size $n$, recent work has explored the incorporation of monotonicity (or isotonicity) information, which can often be elicited from physical systems. In practical applications with many input variables, however, existing Bayesian isotonic models can face statistical and computational limitations, which may result in worse performance compared to models that do not incorporate isotonicity. We propose a new transformed additive isotonic model (TAIM), which aims to tame this "curse-of-dimensionality". TAIM makes use of a flexible transformed additive isotonic modeling framework, which leverages a data-estimated link transformation and a monotone basis model with spike-and-slab priors on basis weights. Prediction-wise, TAIM achieves (up to log factors) a posterior contraction rate of $O(n^{-1/3})$ when the true black-box function is in a transformed additive isotonic form with mild smoothness conditions. Such a rate does not depend on the input dimension $d$ for terms involving $n$, which softens the effect of dimensionality on posterior predictions. Computation-wise, TAIM allows for efficient posterior inference via a carefully designed Gibbs sampler, where each sampling iteration requires only linear work in $d$. We further present an extension of TAIM that can model potential deviations from transformed additivity. Numerical experiments and two applications show the effectiveness of TAIM for isotonic surrogate modeling with many input variables.

Plug-and-Play (PnP) algorithms are a class of iterative methods for inverse imaging. Within an optimization algorithm, they combine a flexible fidelity term, encoding the forward operator, and a pretrained image denoiser, in order to deal with more severe corruptions such as blurring or downsampling when reconstructing an image. This work studies provably convergent PnP methods for mixed-noise forward processes by using the infimal convolution as a fidelity term, providing a statistical interpretation as a joint maximum a-posteriori estimator over both noises, and preserving the Bayesian MAP interpretation of PnP methods. Independently, we extend the PnP formulation to multiple prior terms using the Davis–Yin three-operator splitting. This extension can be combined with either standard fidelities or the proposed mixed-noise infimal-convolution fidelities. We verify that these generalized PnP methods are convergent under standard Kurdyka–Lojasiewicz conditions. Numerical experiments on Laplace-Gaussian and Poisson–Gaussian noise demonstrate stable convergence of single-prior and multiple-prior PnP methods, and divergence under fidelity mismatch. Furthermore, PnP with infimal convolution fidelities are able to scale to noise up to 38% standard deviation, with multiple priors reaching different stationary points that qualitatively preserve more textural properties.

Background. Magnetic resonance imaging (MRI) is central to Alzheimer's disease research and clinical evaluation, yet timing of first MRI after cohort entry varies substantially. In aging cohorts, many participants die before imaging, and treating death as noninformative censoring can distort inference on neuroimaging access and conceal disparities. Objective. To quantify individual-level determinants and center-level heterogeneity in timing of first neuroimaging while accounting for the competing risk of death. Methods. We analyzed National Alzheimer's Coordinating Center Uniform Data Set v3 linked to MRI records from 2015 through March 2025. The study included 20,867 participants aged 50-95 years without dementia at baseline and with known APOE {\epsilon}4 status. Time from baseline to first post-baseline MRI was the primary event, with death before MRI treated as a competing event. We fit Bayesian multilevel cause-specific accelerated failure time models with Alzheimer's Disease Center random intercepts, using Weibull and log-logistic models for MRI timing and death, respectively. Frequentist accelerated failure time and Cox models were used for sensitivity analyses. Results. APOE {\epsilon}4 was not credibly associated with timing of first MRI. Older age, lower educational attainment, and non-White race were associated with delayed MRI acquisition. Center-level variation in MRI timing was substantial. Conversely, APOE {\epsilon}4 homozygosity was associated with earlier death prior to MRI, alongside strong effects of age and cognitive impairment. Conclusions. Timing of first MRI in Alzheimer's research cohorts is driven primarily by demographic, social, and center-level factors rather than genetic risk. Accounting for competing mortality and site heterogeneity is necessary to avoid biased neuroimaging samples and support equitable Alzheimer's research.