Papers

We study the number of exact matrix-vector queries needed to sample a centered $d$-dimensional Gaussian with condition number at most $\kappa$ to total variation error $\varepsilon$. The best known algorithm uses $O(\min\{d,\sqrt{\kappa}\log(\frac{d}{\varepsilon})\})$ queries, while the known lower bounds do not grow as $\varepsilon\to0$ and match this bound only for $\kappa\le d^{1/5-\delta}$ at constant accuracy. We introduce a construction based on a hidden determinant sign and prove a lower bound of $\Omega\!\left(\min\left\{d,\sqrt{\kappa}\log\left(1+\frac{d}{\sqrt{\kappa}\,\varepsilon}\right)\right\}\right)$ uniformly for all $d\ge1$, $\kappa\ge4$, and $0<\varepsilon\le\tfrac{1}{16}$. This bound matches the known upper bound up to constant factors whenever $\kappa\le d^{2-\delta}$ for fixed $\delta\in(0,2)$, or $\varepsilon\le d^{-c}$ for fixed $c>0$. We also extend a recent fixed-spectrum construction with a hidden eigenbasis from $\kappa=2$ to $2\le\kappa\le d^{1-\delta}$, for any fixed $\delta\in(0,1)$. This gives an $\Omega_\delta(\sqrt{\kappa}\log d)$ lower bound at constant accuracy that holds even when the spectrum of the precision matrix is known to the algorithm.

Diffusion models are a powerful framework for learning complex probability distributions. While widely used diffusion models for continuous data rely on Gaussian noising processes, we instead develop diffusion models based on Poisson noise. More precisely, we introduce and study a new mathematical framework, called Mixed Poisson Diffusion Models (MPDM), with noising dynamics driven by multivariate counting processes. Our methodology is motivated by Poisson inverse problems. The success of Gaussian diffusion models has led to their use as prior models in a growing number of posterior sampling methods for inverse problems. However, most of these approaches focus on Gaussian inverse problems, in which observations are generated through a deterministic forward operator corrupted by additive Gaussian noise. By contrast, Poisson inverse problems, where observations are discrete and Poisson distributed, remain comparatively underexplored in this literature. Building on our Poisson generative dynamics, we develop a natural framework for posterior sampling in Poisson inverse problems and derive several posterior sampling algorithms based on the proposed model. We evaluate our approach on realistic Poisson inverse problems and show that the resulting framework is well suited to this class of problems.

It is folklore that passing "convergence" diagnostics does not guarantee that a Markov chain has converged. We prove this as an impossibility for any finite checklist of diagnostics. Running more chains, longer chains, and smarter diagnostics does not avoid this impossibility in certifying convergence. Under regularity assumptions, we determine how the worst-case undetected non-convergence decays as more diagnostics are added, and how to choose diagnostics optimally. Folklore suggests that the sampler's slowest mixing modes are the "dangerous" ones that must be monitored; we prove that the slow modes yield nearly optimal diagnostics as sampling proceeds. Given a source class, we also give an anytime-valid bound on the remaining discrepancy and show how undetected bias caps the effective sample size that a finite checklist can certify.

The unadjusted Langevin algorithm (ULA) and its stochastic-gradient variant, SGLD, are workhorses of modern Bayesian machine learning, but can be fragile when faced with a target distribution possessing a non-Lipschitz drift. In such settings, the expected squared jump distance can grow rapidly, leading to numerical instability or the need for prohibitively small step-sizes. We propose the saturated unadjusted Langevin algorithm (SULA) that pairs a saturated gradient drift with a gradient-informed state-dependent diffusion, while preserving first-order consistency with the Langevin diffusion. The drift and diffusion terms are chosen to guarantee that the one-step squared update is independent of the gradient. Consequently, SULA is robust to the choice of step-size and preserves estimation and prediction quality. We provide sufficient conditions under which SULA with a step-size $h$ converges weakly to the Langevin diffusion at the standard $O(h)$ rate, the SULA chain is geometrically ergodic, and that the bias of its invariant measure relative to the target is $O(h)$. Further, we develop SG-SULA, a stochastic-gradient version of SULA that rescales the stochastic gradient using a Gaussian CDF correction under Gaussianity of the gradient noise. Numerical experiments on various examples demonstrate SULA's robustness to tuning and SG-SULA's improved predictive performance.

Monte Carlo estimation, an essential tool for quantifying uncertainty in complex scientific and engineering applications, often has unacceptably large errors when computing resources are limited and each sample requires evaluating an expensive high-fidelity model. Multifidelity control variate methods seek to increase the accuracy of Monte Carlo estimation by combining a limited number of expensive high-fidelity samples with a larger number of cheaper lower-fidelity samples. Many different multifidelity estimators have been proposed, including multilevel Monte Carlo (MLMC), multifidelity Monte Carlo (MFMC), approximate control variates (ACVs), and multilevel best linear unbiased estimators (MLBLUEs). An important design choice in multifidelity estimation is the allocation of a limited computational budget among the available models to minimize estimator error. The state of the art in solving this budget allocation problem relies on numerical optimization; to date, closed form solutions have only been derived for estimator forms known to be suboptimal. In this work, we provide closed form expressions for optimal fractional sample allocations for a large class of nested multifidelity control variate estimators, in which higher-fidelity models are evaluated on a subset of the inputs that lower-fidelity models are evaluated on. Our analysis is based on reformulating nested MLBLUEs as an extended version of the MFMC estimator, and is applicable to general multivariate Monte Carlo estimation. We additionally propose a new approach for integer-valued budget allocation based on a greedy strategy used in combinatorial optimization, as well as a method for robust approximation of control variate coefficients from limited pilot samples. Two numerical examples illustrate our theoretical results and demonstrate the increased accuracy of our extended MFMC formulation relative to existing estimators.

Functional brain aging is often characterized by differences between people of different ages, although aging unfolds within individuals. This distinction is especially important for functional connectivity, where the direction of covariance change may evolve across adulthood. We introduce a Bayesian longitudinal covariance-operator framework that models each visit covariance as a latent positive-definite state informed by repeated resting-state fMRI runs. The approach calibrates run-level information for temporal dependence and replicate variability, separates within-person aging from between-person age differences, and uses an adaptive tangent representation to learn the low-dimensional structure of longitudinal change without relying on pooled covariance variation. The longitudinal operator varies smoothly in magnitude and orientation with age, while exact Fréchet derivatives translate the fitted trajectory to the covariance scale. This framework allows us to examine how brain covariance organization changes within individuals, whether longitudinal aging follows the population age gradient, and which functional systems contribute to their divergence. In the Dallas Lifespan Brain Study, within-person covariance aging differs substantially from the cross-sectional age gradient through early and middle adulthood, with the discrepancy diminishing later in life. Covariance reorganization is strongest around midlife, and the longitudinal signal is concentrated in a few coordinated directions, with most off-diagonal change occurring between functional systems. These findings demonstrate that population-level signatures of brain aging do not necessarily reflect individual aging trajectories and highlight the importance of modeling longitudinal covariance change directly.

Reliable uncertainty quantification remains an important challenge in AI. Hierarchical Bayesian methods provide a principled framework for this by representing unknown quantities through posterior distributions. Joint inference of neural-network parameters and prior hyperparameters is challenging for conventional Markov chain Monte Carlo (MCMC) methods due to strong dependencies. In this work, we introduce the multilevel particle marginal Metropolis-Hastings (MLPMMH) method, which combines pseudo-marginal inference for the hyperparameters with coupled sequential Monte Carlo (SMC) across hierarchy of network-width truncations for conditional parameter inference. The hierarchy is induced by a trace-class neural network (TNN) prior with width-dependent decay, motivated empirically by the compressible structure observed in trained neural-network parameters. Under standard multilevel coupling assumptions, we establish that MLPMMH achieves mean-square error (MSE) $\mathcal{O}(\varepsilon^2)$ at expected computational cost $\mathcal{O}(\varepsilon^{-2})$, and verify the predicted rate on a Bayesian regression problem. To address larger-scale applications, we further introduce S-MLPMMH, a practical approximation that replaces the exact sampling components with a PCA-based hierarchy and scalable Bayesian Monte Carlo. Experiments on MNIST demonstrate competitive predictive performance and uncertainty quantification while exploiting the parallel structure of the approximation. Software has been made publicly available at https://github.com/submission-code-anon/mlpmmh-anon.

We consider the problem of estimating the expectation of a function (observable) with respect to a high-dimensional probability distribution, under the assumption that a low-dimensional collective variable map captures all the directions of slow mixing for the distribution. More precisely, under the assumption that the fibers of the collective variable map satisfy a log-Sobolev inequality, we establish non-asymptotic relative-error guarantees for this problem using the eigenvector method for umbrella sampling (EMUS). Our analysis applies to both log-Lipschitz observables and indicator of convex sets, and gives complexity bounds that are polynomial in regularity parameters and inverse relative accuracy, and exponential in the dimension of the collective variable space. As a relative error guarantee, our bounds apply to rare event estimation. We consider Gaussian biasing windows defined on a low-dimensional collective variable space and analyze a modified sparse matrix version of EMUS. The bounds quantify the effect of the map through mixing of the conditional distribution and fluctuation of the observable along the fibers.

We construct an approximate sampler for finite doubly free-boundary Macdonald rail-yard processes with Macdonald, boundary, and edge parameters in $(0,1)$. The construction allows arbitrary finite sequences of additions and removals of horizontal or vertical strips and all four boundary types. It reverses a finite commutation–reflection history from a partition sequence whose entries are equal, with the common partition weighted by both boundaries. We derive local reconstruction kernels from the skew Cauchy, Littlewood, and branching identities. Bulk and boundary updates combine integer increments with conditional shape sampling at a prescribed size; same-sign swaps use finite conditional resampling. For $K$ reflection periods and an initial partition-size cutoff $M$, the algorithm terminates almost surely and samples the resulting approximating law exactly. We prove total variation convergence to the target process as $K,M\to\infty$ and give explicit bounds for both truncation errors, decaying geometrically in the cutoffs for fixed input parameters. Under fixed parameter bounds, we also bound the expected operation count of direct enumeration in a unit-cost model. We further construct an exact rejection sampler for the initial partition law without a size cutoff or evaluation of its normalizing constant. Using this sampler removes the size-truncation error, leaving only the error from the finite reflection history.

Sampling from high-dimensional, heavy-tailed distributions is a fundamental challenge in computational statistics, as many standard Markov chain Monte Carlo (MCMC) methods mix poorly in such settings. Recently, Stereographic MCMC [Yang et al., 2024] and the Sub-Cauchy Projection Sampler [Grazzi et al., 2026] have been shown to perform successfully on such tasks. However, establishing their non-asymptotic convergence properties remains an important open problem. In this work, we fill this gap by establishing non-asymptotic upper bounds on the mixing time of the stereographic projection and sub-Cauchy projection samplers. Our results demonstrate that, under certain conditions on the target and initial distributions, the mixing time is polynomial in dimension for a broad class of distributions, including light- and heavy-tailed cases. Motivated by the theoretical analysis, we further establish a new weighted isoperimetric inequality that extends the classical version for log-concave distributions to the heavy-tailed setting with optimal dimension dependence. The proof techniques provide new insights into the geometric properties of heavy-tailed distributions that govern rapid mixing in high dimensions.

Coordinate descent solvers for regularized generalized linear models (GLMs) differ in how they update the intercept. This seemingly routine choice can substantially affect convergence when the response is imbalanced. For direct coordinate descent on the original GLM loss, we compare a conservative step based on worst-case curvature, a safeguarded Newton step based on current curvature, and repeated Newton steps that solve the intercept subproblem. The conservative update can leave much of the interaction between the intercept and the coefficients unresolved, producing a slowdown when that interaction controls convergence. A single Newton step avoids this attenuation and, in our experiments, follows nearly the same outer trajectory as exact minimization while requiring fewer inner steps. Simulated and real-data experiments, including comparisons across six production solvers and a controlled intervention within skglm, support this explanation. We also show how the same mechanism appears in quadratic-approximation methods and extends to multinomial models with rare classes. We recommend a safeguarded Newton update for direct coordinate descent and full intercept optimization within frozen local quadratic surrogates.

Sampling Gaussian measures restricted to nonlinear manifolds in high dimensions requires a proposal kernel that respects both geometric constraints and dimensional scaling requirements. We introduce tangent preconditioned Crank-Nicolson (pCN), which combines a manifold projection with a cancellation of the ambient Gaussian contribution, leaving only a geometric residual. Considering a shifted sphere with $\sqrt d$ scaling of radius and centre as a tractable benchmark embedded manifold, we prove that uncorrected tangent pCN results in a proposal noise scale of order $d^{-1/6}$ (compared to $d^{-1/2}$ for a tangent random-walk). We establish asymptotic optimal tuning for the geodesic expected squared jump distance which is maximised at an acceptance probability of 0.574 for the family of drift-corrected tangent proposals containing tangent pCN and tangent MALA. We also prove a joint diffusion limit for the parallel and orthogonal coordinates to the sphere's centre vector of tangent pCN. The orthogonal diffusion speed matches the limiting geodesic expected squared jump distance and is maximised at acceptance probability 0.574, whereas the parallel speed is maximised at 0.712, demonstrating that optimal tuning depends on the observable. We then identify an intrinsic precision correction characterised by potential Hessian and Ricci curvature that yields a proposal kernel with dimension-independent fixed step size limit. We discuss a generalisation to arbitrary hypersurfaces and present nondegenerate scaling of tangent pCN applied to conditioned diffusions under resolution refinement.

In modern sequential Monte Carlo (SMC) applications, such as diffusion model samplers, propagating a particle is expensive, so only a few particles can be afforded. In this regime, each resampling step decides which regions of the state space survive, so a better resampling scheme is worth extra computation. However, classical resampling methods only replicate existing particles and suffer from sample impoverishment, especially when the particle size is small relative to the dimension. To create new particles, we relax the unbiasedness of resampling and control the mean squared error via a Wasserstein criterion. Specifically, we propose Wasserstein-Dirichlet resampling (WDR), which constructs an equal-weight measure that optimally approximates the weighted particle distribution by solving a free-support Wasserstein barycenter problem. We then randomize the resulting optimal coupling by Dirichlet mixing. Theoretically, we provide non-asymptotic mean squared error bounds that quantify the bias-variance trade-off and show that the rate depends on the intrinsic dimension. A truncated variant preserves high-weight particles while applying WDR to the remaining mass, targeting mild weight degeneracy. Experiments on state-space models and conditional image generation demonstrate that WDR consistently improves estimation accuracy, stability, and sample diversity. The improvements persist at matched wall-clock cost when the propagation cost dominates.

The Hybrid forecast covariance model is an operationally successful way to build a background error covariance matrix. It is constructed as a weighted sum, with scalar weights, of the climatological covariance matrix and a localized ensemble covariance matrix, with a localization matrix typically depending on a localization length. For best performance, the Hybrid needs to be optimized by tuning those weights and, possibly, other parameters, such as the localization length. Here, a systematic and automated method for this tuning operation is introduced. This technique is referred to as hypervariational method, and it is formulated in terms of a cost function, built using a time series of innovations as an essential ingredient. This cost function, which depends on the parameters, is minimized when they have optimal values. Equations for its minimization, based on Newton's method, are derived. Besides the aforementioned, commonly used Hybrid, the hypervariational method is applicable to any known Hybrid (and to any set of parameters on which it depends analytically), including recently developed Hybrids, not necessarily constructed as weighted sums. Here, the hypervariational method is applied to tune the ensemble covariance matrix weight in the seminal-formulation Hybrid (which does not feature localization). This is done in a synthetic experimental setting where a theoretically-derived expression for the optimal value for that weight is known. Moreover, in that experimental setting, a further method, based on a linear fit procedure, to determine that optimal value is possible. Both the hypervariational method and that linear fit procedure are shown to produce an optimal weight value matching (up to a tolerance of $0.1$) the optimal value known from the theory. Furthermore, Newton's minimization in the hypervariational method is shown to converge immediately (at the first recursive step).

Bridge sampling estimates of marginal likelihoods are routinely reported with a delta-method Monte Carlo standard error (MCSE). Existing theory is asymptotic and cannot say whether a given run is long enough to trust the estimate or its MCSE. We give finite-sample answers. The optimal bridge estimator fits a logistic regression that separates posterior from proposal draws, and its accuracy is governed by their overlap, within a factor of two of the classification's Bayes error. Because the estimating function is a sum of bounded terms, no moment conditions are needed, and all bounds hold uniformly over targets, proposals and dimensions. A number of draws inversely proportional to the overlap suffices and is necessary in the worst case, paralleling the exponential-in-Kullback–Leibler requirement of importance sampling. A Berry–Esseen bound with an explicit constant holds at the optimal rate, and the MCSE is accurate once the effective overlap is large; it is estimated at no cost. For Markov chain draws, the required sample size grows in proportion to the sampler's relaxation time, unavoidably in the worst case. The estimated bridge, the fixed-point iteration and the fitted proposal matter only at second order, and we give free diagnostics for sampler failure and misestimated autocorrelation.

We present CLASSO, a path-following method for the LASSO subject to general linear equality and inequality constraints. The local direction is obtained from a quadratic programme on an augmented active set, which determines both the coefficients that leave zero and the weakly active inequalities that remain binding. For a positive-definite Hessian and a nonempty feasible set, this programme has a unique primal minimiser, even when the binding constraint rows are dependent, and gives the local derivative of the solution path. The construction recovers the LASSO-modified LARS direction without constraints, and over the full cone gives the segment leaving the origin. We prove that the primal path is piecewise affine with finitely many kinks, without a constraint qualification or unique multipliers. A finite path does not imply a finite traversal: re-certifying at each kink can fail to terminate. CLASSO instead follows the minimum-norm certificate along a lexicographic rate; we prove that this rule reaches $\lambda=0$ in exact arithmetic within an explicit bound and traverses the whole path. These results are machine-checked in Lean~4, and the floating-point implementation agrees with an exact rational one on all $210$ small test problems. We assess the implementation on a $39$-case synthetic suite and six constrained versions of the diabetes benchmark. On the $31$ nontrivial synthetic cases with $n\ge p$, sampled relative objective differences from numerical reference solves are below $10^{-9}$ in magnitude; unnormalised endpoint differences reach $1.5\times10^{-7}$ in magnitude. For $n<p$, matched differences range from $-3.6\times10^{-7}$ to $2\times10^{-9}$ across five ridge strengths. Across the synthetic and diabetes benchmarks, checks of all $580$ recorded segments found no missed interior event and all checked KKT residuals within the stated tolerances.

This article develops optimal order (up to logarithmic factors) weighted cubature rules for Monte Carlo numerical integration in Sobolev reproducing kernel Hilbert spaces over common matrix manifolds. Specifically, we consider matrix models of the Stiefel, Grassmann, flag, and Lagrangian Grassmann manifolds each over either real numbers $\mathbb{R}$, complex numbers $\mathbb{C}$, or quaternions $\mathbb{H}$. The proposed cubature rules use evaluation points which are independent and identically distributed uniform on the manifold. The weights are chosen to be optimal with respect to the reproducing kernel Hilbert space worst-case error, under the constraint that the weights sum to one. Crucially, this constraint ensures optimal weights may be found without having to evaluate the kernel mean embedding (the integral of the kernel), which is often difficult to compute. We consider computable Matérn kernels of Frobenius distances between matrices which reproduce the Sobolev reproducing kernel Hilbert spaces with equivalent norms. Open-source numerical experiments validate our theoretical development, including an application to simultaneous approximation of intrinsic volumes in the Cauchy–Kubota formula. We also present preliminary empirical results beyond independent points including transformations of low-discrepancy quasi-Monte Carlo points from the unit cube and jointly optimized points and weights.

The multivariate linear mixed model combined with restricted maximum likelihood estimation of (co)variance components and solving the mixed model equations is widely used for the analysis of plant breeding data. Several algorithms for restricted maximum likelihood estimation have been proposed, with the average information algorithm being widely used for its computational efficiency, quadratic rate of convergence, and ability to handle large amounts of parameters. A key reason for the scalability of linear mixed models is the exploitation of sparse matrix techniques. Recently, decreasing costs of high-throughput data collection have led to high-dimensional sets of features characterizing genotypes and environments. A convenient way of using these features is through the use of kernel-based covariance structures. However, while convenient, these kernels typically have dense inverses meaning optimizations based on sparse matrix techniques are not always applicable. For these settings, the introduction of general-purpose computing on graphics processing units brought significant computational advantages for dense linear algebra operations. We present a hybrid GPU/CPU implementation of the average information restricted maximum likelihood algorithm. We eliminate sparse equations using supernodal Cholesky factorization and rely on CUDA GPU acceleration to solve the resulting reduced dense system. We show the significant computational advantages of this approach for fully dense problems, as well as problems with a mix of dense and sparse equations.

Exact reversible locally adaptive step sizes are introduced by Bou-Rabee et al. (2026) in the context of the No-U-Turn Hamiltonian Monte Carlo (HMC) sampler for the purpose of high performance HMC sampling in challenging target distributions. This paper carefully explains and illustrates the locally adaptive step size procedure under more general assumptions. Alternative error measures are introduced with the aim of further enhancing MCMC performance by reducing the prevalence of reversibility failure rejections. In addition, the paper shows how to introduce the locally adaptive step sizes also within standard HMC and Metropolis-corrected kinetic Langevin/generalized HMC (GHMC). Finally, a numerical study of the relative merits of the proposed error measures in the context of different Hamiltonian flow-based MCMC methods is provided.

Existing cumulative regret guarantees for Bayesian optimization (BO) do not necessarily characterize the quality of the algorithm's final output. This gap is critical under noisy evaluations, where the goal is to return a single high-quality recommendation. To address it, we introduce efficient exploration strategies that can be incorporated into any BO procedure by occasionally replacing the acquisition query with an exploration query. To validate this approach, we bound the Gaussian process posterior standard deviation over the search space via the fill distance of the queries, which measures the largest region left unexplored. We then use this result to derive upper bounds on the inference regret, i.e., the gap between the global optimum and the function value at the algorithm's returned point. Empirically, these exploration strategies improve the final recommendations of existing BO algorithms on synthetic benchmarks, hyperparameter optimization, and a control task.

Split Gibbs sampling (SGS) is a versatile tool to sample from complex posterior distributions, such as those arising in imaging inverse problems. However, existing SGS instances are limited to quadratic variable couplings, restricting their applicability to Gaussian likelihoods. This letter genearalizes SGS to non-Euclidean geometry by introducing a coupling kernel defined through a Jensen-Bregman divergence. This framework not only enables posterior sampling in these settings but also unifies several ad hoc sampling schemes from the literature under a principled theoretical foundation.

Verifier-guided inference-time scaling uses additional computation to score intermediate states and steer generation toward better outputs. However, imperfect intermediate verifiers can distort sampling, and repeated refinement alone does not ensure recovery of the desired reward-tilted distribution. We introduce Continuous-State Verifier-Guided Backtracking (CS-VGB), an inference-time scaling method for sampling from reward-tilted distributions with continuous generative models through repeated denoising and renoising. With correctly paired denoising and renoising transitions, the stationary distribution of CS-VGB at the clean checkpoint equals the desired reward-tilted target, even when intermediate verifiers are inaccurate. Under stated regularity conditions, CS-VGB converges exponentially to the target despite verifier error, and the expected number of updates needed to reach a clean sample grows only linearly with the number of checkpoints. In analytically tractable settings, we empirically demonstrate that CS-VGB recovers the desired reward-tilted distribution despite verifier error, whereas the compared verifier-based sampling methods fail to recover the target under the same verifier errors. We also instantiate the framework with diffusion models and one-step stochastic posterior maps. Across prompt-guided ImageNet generation, facial-attribute generation on CelebA, and multi-robot planning, recurrent CS-VGB search improves reward and task success as inference budgets grow, with image-generation rewards continuing beyond the observed plateaus of denoising-only baselines.

Generalized neural network mixed models (GNMMs) combine neural networks and mixed-effects models to analyze longitudinal and clustered data. Although GNMMs have shown good predictive performance in practice, their statistical properties are not well understood. Therefore, this paper studies several important properties of GNMM and provides theoretical support for its use. First, we show that the neural-network component can approximate an unknown nonlinear fixed-effect function arbitrarily well and establish consistency of the fitted neural- network function for its population target within a given network class. Next, we establish consistency of random-effect estimation under increasing within-subject information and show that the GNMM predictor converges to the true predictor. We then prove consistency of the conditional mean function and show that the excess mean squared prediction error (MSPE) relative to the oracle prediction error converges to zero. We conduct simulation studies to examine nonlinear fixed-effect recovery, random-effect estimation, predictor and conditional mean estimation, and excess prediction error under different numbers of subjects and repeated observations. We also apply GNMM to a longitudinal schizophrenia dataset to predict continuous and binary depression outcomes using smartphone- based behavioral measurements and compare its predictive performance with generalized linear mixed models and conventional artificial neural networks. Both the theoretical and empirical results demonstrate that GNMM provides a flexible framework for modeling nonlinear longitudinal and clustered data. These findings provide a statistical foundation for GNMM prediction and support future methodological developments.

Clustering is typically performed using only object features, although additional information about relationships between objects may be available and can be incorporated through constraints. Among these, pairwise constraints require two objects to belong to the same cluster or to different clusters. Several extensions of the k-means algorithm accommodate such constraints, but they sometimes fail to find a feasible assignment or incur high computational cost. Motivated by these limitations, we propose an algebraic approach to k-means clustering with pairwise constraints. The algorithm simplifies the constraints using Gröbner bases before optimization, improving the efficiency of cluster assignment. It can also detect whether the constraint set has a feasible assignment before optimization. Comparative experiments with existing algorithms showed that the proposed algorithm achieved comparable or better clustering quality at a lower computational cost, particularly when the numbers of constraints and clusters were small.

Human reference labels used to evaluate medical AI can vary across readers and annotation workflows. We quantified agreement of a fixed pathology AI system with a future qualified reader without treating inferred consensus as error-free biological truth. Linked Evaluation and Unbiased NuCLS multirater data were analyzed using a fixed external HoVer-Net decision vector. The selected frame contained 1,144 anchors, 52 fields of view, five slides/patients, and five primary readers. The target was fixed-AI agreement with a future reader, with descriptive contrast $\Delta=\theta_E-\theta_U$. Analyses included raw agreement and iMRMC benchmarks, reader-only and shared-detectability Bayesian models, structural sensitivities, adversarial simulations, and predictive assessment. Among 5,188 reader-anchor pairs observed in both conditions, 19.62% changed AI-agreement state, while raw agreement was 0.5466 in Unbiased and 0.5473 in Evaluation. The reader-only model gave $\Delta=-0.0008$ (95% CrI, $-0.0033$ to $0.0008$), and the shared-detectability model gave $\Delta\approx0$ ($-0.0140$ to $0.0127$). Across structural perturbations, absolute agreement ranged from 0.420 to 0.646 in Unbiased and 0.409 to 0.640 in Evaluation. Frozen adversarial simulations met the comparative-stability criterion (median $R=0.295$; 14/16 model-scenario cells had $R<1$). Held-out-AI prediction supported shared predictive information, but varied by slide composition. The Evaluation-minus-Unbiased contrast was comparatively stable, whereas absolute agreement was structurally sensitive. These findings do not establish equivalence, causal assistance effects, biological accuracy, or new-patient transportability. Evaluations using variable human references should report reader population, selected frame, observation mechanism, and structural sensitivity.

The Copula Active Subspaces (CAS) method of Part I estimates a non-Gaussian observation-noise density $\pi_{\boldsymbol{n}}$ from residual samples. Its KL divergence $\mathcal{D}_{\mathrm{KL}}(\pi_{\boldsymbol{n}}\,\|\,\widehat\pi_{\boldsymbol{n}})$ admits an exact decomposition into marginal, projection and reduced-density terms and bounds the posterior KL divergence averaged over the true predictive distribution. We derive a priori bounds for these terms, together with estimates or bounds from residual samples. Weighted nearest-neighbour entropy estimates are used for the marginal and reduced-density terms, with asymptotic standard errors in exact rank coordinates. For the projection error, the score is fitted on one subsample, and its mean squared magnitude outside the retained subspace is estimated on the other. Corrections for score and rank-transform errors then give a theoretical upper bound with probability at least $1-\delta$. Tests on three known noise laws in dimension $d=20$ use the estimated rank transform. In $11$ of $12$ combinations of example and sample size, the sum of the two estimated terms and the trace bound exceeds the reference divergence after averaging over the estimated models. In these tests, score-estimation error is the largest contribution to looseness. The overestimate from the Gaussian log-Sobolev inequality alone is at most a factor of $1.98$.

The cvxoptdes R-package provides a unified framework to perform convex optimal experiment design and online response optimization for linear, generalized linear and additive, nonlinear, and random intercept mixed models, using exclusively open-source solvers. It supports optimization and augmentation of weighted, penalized and convex constrained continuous designs across a discrete grid of candidate points and rounding to integer-valued designs using point-exchange algorithms or mixed-integer programming for orthogonally constrained designs, including the class of orthogonal minimally aliased response surface (OMARS) designs. Models and design constraints are defined and updated in an object-oriented context with the aim to facilitate the process of sequential (penalized) design augmentation analogous to online design procedures known from active learning and Bayesian optimization.

We present a simulation-based method for the estimation, sensitivity analysis, and optimization of distortion risk measures (DRMs), a class of performance measures that capture a wide range of risk attitudes through an underlying distortion function. Existing simulation approaches often approximate DRMs using numerical partitions or empirical distributions. We show that, when the distortion function is piecewise polynomial, the nonlinear structure of a DRM can be handled exactly, without introducing distortion discretization error, using only simulated system outputs. This representation enables the development of a computational method that integrates finite-difference techniques for sensitivity analysis and for solving a broad class of differentiable risk-based optimization problems. Under suitable conditions, we establish strong consistency and convergence rates for DRM estimation, sensitivity analysis, and optimization. We further extend the approach to non-polynomial distortion functions through polynomial approximation. Numerical experiments on synthetic and calibrated financial examples demonstrate the effectiveness and robustness of the proposed approach across a range of applications.

Sensitivity analysis asks how strong unmeasured confounding would have to be before an observational conclusion changes. Established methods answer this question using different quantities, including partial $R^2$ values, correlation products, and risk ratios. In this article, we present confoundvis, an R package that organizes these results through a common object, plotting functions, and reports. It calculates partial $R^2$ and impact-threshold paths from fitted linear models, imports results from sensemakr, konfound, and EValue, and provides measured-covariate benchmarks for the OLS frameworks. We demonstrate the workflow with 5,735 patients from the right heart catheterization study. The adjusted linear probability coefficient is 0.0634, the point-null robustness value is 0.0608, and the conditional significance impact threshold is 0.0376. The largest adverse leave-one-term-out measured-covariate product is 0.0090. Because the benchmark and threshold use different conditioning sets, their ratio is descriptive. A separate modified Poisson model gives a risk ratio of 1.103 and an E-value of 1.439. The contribution is a software workflow that keeps each method's units and decision rule clear. The package is available from CRAN.

Latent field models are essential for spatial statistics, but the posterior distribution of the latent field, even with a Gaussian process (GP) prior, is analytically intractable whenever the likelihood is non-Gaussian. Existing approximations are either too slow or too restrictive for large non-Gaussian datasets. We present a scalable variational inference framework that mirrors the sparse nearest-neighbor autoregressive structure of the prior in the variational family, uses normalizing-flow conditionals for expressiveness, and amortizes all variational parameters via a single shared neural conditioner. Doubly stochastic optimization with reduced ancestor sets achieves polylogarithmic complexity per iteration.