Papers

Diffusion-based solvers for partial differential equations (PDEs) are often bottle-necked by slow gradient-based test-time optimization routines that use PDE residuals for loss guidance. They additionally suffer from optimization instabilities and are unable to dynamically adapt their inference scheme in the presence of noisy PDE residuals. To address these limitations, we introduce PRISMA (PDE Residual Informed Spectral Modulation with Attention), a conditional diffusion neural operator that embeds PDE residuals directly into the model's architecture via an attention-inspired modulation mechanism in the spectral domain, enabling gradient-descent free inference. In contrast to previous methods that use PDE loss solely as external optimization targets, PRISMA integrates PDE residuals as integral architectural features, making it inherently fast, robust, accurate, and free from sensitive hyperparameter tuning. We show that PRISMA has competitive accuracy, at substantially lower inference costs, compared to previous methods across five benchmark PDEs especially with noisy observations, while using 10x to 100x fewer denoising steps, leading to 15x to 250x faster inference.

Motivated by an application to study the impact of temperature, precipitation and irrigation on soybean yield, this article proposes a sparse semi-parametric functional quantile model. The model is called “sparse” because the functional coefficients are only nonzero in the local time region where the functional covariates have significant effects on the response under different quantile levels. To tackle the computational and theoretical challenges in optimizing the quantile loss function added with a concave penalty, we develop a novel Convolution-smoothing based Locally Sparse Estimation (CLoSE) method, to do three tasks in one step, including selecting significant functional covariates, identifying the nonzero region of functional coefficients to enhance the interpretability of the model and estimating the functional coefficients. We establish the functional oracle properties and simultaneous confidence bands for the estimated functional coefficients, along with the asymptotic normality for the estimated parameters. In addition, because it is difficult to estimate the conditional density function given the scalar and functional covariates, we propose the split wild bootstrap method to construct the confidence interval of the estimated parameters and simultaneous confidence band for the functional coefficients. We also establish the consistency of the split wild bootstrap method. The finite sample performance of the proposed CLoSE method is assessed with simulation studies. The proposed model and estimation procedure are also illustrated by identifying the active time regions when the daily temperature influences the soybean yield.

Modeling heterogeneous correlated time series requires the ability to learn hidden dynamic relationships between component time series with possibly varying periodicities and generative processes. To address this challenge, we formulate and evaluate a windowed variance-correlation metric (WVC) designed to quantify time-varying correlations between signals. This method directly recovers hidden relationships in an specified time interval as a weighted adjacency matrix, consequently inferring hidden dynamic graph structure. On simulated data, our method captures correlations that other methods miss. The proposed method expands the ability to learn dynamic graph structure between significantly different signals within a single cohesive dynamical graph model.

An important class of spatio-temporal models is constructed by leveraging the hierarchical structure of dynamical (or, state-space) models. This paper proposes a new statistical dynamical model for spatio-temporal processes motivated by second-order stochastic partial differential equations (SPDE). In particular, an infinite-dimensional linear state-space representation is obtained where the state transition is governed by a proposed SDE. Then, using the Galerkin's method, a finite-dimensional approximation to the infinite-dimensional SDE is obtained, yielding a dynamical model with finite states that facilitates computation and parameter estimation. The space-time covariance of the approximated dynamical model is obtained, and the error between the approximate and exact covariance matrices is quantified. Comprehensive numerical investigations, including 2D wave equation, seismic wave propagation, advection-diffusion equations and wildfire aerosol propagation processes, are performed to demonstrate the application of the proposed model. Code is available.

We consider change point detection for the volatility in second order linear parabolic stochastic partial differential equations based on high frequency spatio-temporal data. We give a test statistic to detect changes in the volatility based on change point analysis for diffusion processes and derive the asymptotic null distribution of the test statistic. We also show that the test is consistent. Moreover, we provide some examples and then perform numerical simulations of the proposed test statistic.

Online learning is the cornerstone of applications like recommendation and advertising systems, where models continuously adapt to shifting data distributions. Model training for such systems is remarkably expensive, a cost that multiplies during hyperparameter search. We introduce a two-stage paradigm to reduce this cost: (1) efficiently identifying the most promising configurations, and then (2) training only these selected candidates to their full potential. Our core insight is that focusing on accurate identification in the first stage, rather than achieving peak performance, allows for aggressive cost-saving measures. We develop novel data reduction and prediction strategies that specifically overcome the challenges of sequential, non-stationary data not addressed by conventional hyperparameter optimization. We validate our framework's effectiveness through a dual evaluation: first on the Criteo 1TB dataset, the largest suitable public benchmark, and second on an industrial advertising system operating at a scale two orders of magnitude larger. Our methods reduce the total hyperparameter search cost by up to 10$\times$ on the public benchmark and deliver significant, validated efficiency gains in the industrial setting.

We present Randomized-Accelerated FEAST (RA-FEAST), a hybrid algorithm that combines contour-integration-based eigensolvers with randomized numerical linear algebra techniques for efficiently computing partial eigendecompositions of large-scale matrices arising in statistical applications. By incorporating randomized subspace initialization to enable aggressive quadrature reduction and truncated refinement iterations, our method achieves significant computational speedups (up to 38x on sparse graph Laplacian benchmarks at n = 8000) while maintaining high-accuracy approximations to the target eigenspace. We provide a probabilistic error bound for the randomized warmstart, a stability result for inexact FEAST iterations under general perturbations, and a simple complexity model characterizing the trade-off between initialization cost and solver speedup. Empirically, we demonstrate that RA-FEAST can be more than an order of magnitude faster than standard FEAST while preserving accuracy on sparse Laplacian problems representative of modern spectral methods in statistics.

Many online learning algorithms, including classical online PCA methods, enforce explicit normalization steps that discard the evolving norm of the parameter vector. We show that this norm can in fact encode meaningful information about the underlying statistical structure of the problem, and that exploiting this information leads to improved learning behavior. Motivated by this principle, we introduce Implicitly Normalized Online PCA (INO-PCA), an online PCA algorithm that removes the unit-norm constraint and instead allows the parameter norm to evolve dynamically through a simple regularized update. We prove that in the high-dimensional limit the joint empirical distribution of the estimate and the true component converges to a deterministic measure-valued process governed by a nonlinear PDE. This analysis reveals that the parameter norm obeys a closed-form ODE coupled with the cosine similarity, forming an internal state variable that regulates learning rate, stability, and sensitivity to signal-to-noise ratio (SNR). The resulting dynamics uncover a three-way relationship between the norm, SNR, and optimal step size, and expose a sharp phase transition in steady-state performance. Both theoretically and experimentally, we show that INO-PCA consistently outperforms Oja's algorithm and adapts rapidly in non-stationary environments. Overall, our results demonstrate that relaxing norm constraints can be a principled and effective way to encode and exploit problem-relevant information in online learning algorithms.

Mean-field games (MFGs) study the Nash equilibrium of systems with a continuum of interacting agents, which can be formulated as the fixed-point of optimal control problems. They provide a unified framework for a variety of applications, including optimal transport (OT) and generative models. Despite their broad applicability, solving high-dimensional MFGs remains a significant challenge due to fundamental computational and analytical obstacles. In this work, we propose a particle-based deep Flow Matching (FM) method to tackle high-dimensional MFG computation. In each iteration of our proximal fixed-point scheme, particles are updated using first-order information, and a flow neural network is trained to match the velocity of the sample trajectories in a simulation-free manner. Theoretically, in the optimal control setting, we prove that our scheme converges to a stationary point sublinearly, and upgrade to linear (exponential) convergence under additional convexity assumptions. Our proof uses FM to induce an Eulerian coordinate (density-based) from a Lagrangian one (particle-based), and this also leads to certain equivalence results between the two formulations for MFGs when the Eulerian solution is sufficiently regular. Our method demonstrates promising performance on non-potential MFGs and high-dimensional OT problems cast as MFGs through a relaxed terminal-cost formulation.

Gaussian-process state-space models (GP-SSMs) provide a flexible nonparametric alternative for modeling time-series dynamics that are nonlinear or difficult to specify parametrically. While the Kalman filter is effective for linear-Gaussian trend and seasonal components, many real-world systems require more expressive representations. GP-SSMs address this need by learning transition functions directly from data, while particle filtering enables Bayesian state estimation even when posterior distributions deviate from Gaussianity. This paper develops a particle-filtering framework for GP-SSM inference and compares its performance with the Kalman filter in trend extraction and seasonal adjustment. We further evaluate nonlinear signal-extraction tasks, demonstrating that GP-SSMs can recover latent states under sharp or asymmetric dynamics. The results highlight the utility of combining GP modeling with sequential Monte Carlo methods for complex time-series analysis.

Clinical study populations often differ meaningfully from the broader populations to which results are intended to generalize. Weighting methods such as inverse probability of sampling weights (IPSW) reweight study participants to resemble a target population, but the accuracy of these estimates depends heavily on how well the chosen population represents the population of substantive interest. We conduct a simulation study grounded in empirical covariate distributions from several real-world data sources spanning a continuum from highly selective to broadly inclusive populations. Using treatment effect scenarios with varying levels of effect modification, we evaluate IPSW estimators of the population average treatment effect (PATE) across multiple candidate target populations. We quantify the bias that arises when the dataset used to operationalize the target population differs from the intended inference population, even when IPSW is correctly specified. Our results show that bias increases systematically as target populations diverge from a well-representative population, and that weighting to a poorly aligned target can introduce more bias than not weighting at all. These findings highlight that selecting an appropriate target population dataset is a critical design choice for valid generalization.

Foundation models, and in particular large language models, can generate highly informative responses, prompting growing interest in using these ”synthetic” outputs as data in empirical research and decision-making. This paper introduces the idea of a foundation prior, which shows that model-generated outputs are not as real observations, but draws from the foundation prior induced prior predictive distribution. As such synthetic data reflects both the model's learned patterns and the user's subjective priors, expectations, and biases. We model the subjectivity of the generative process by making explicit the dependence of synthetic outputs on the user's anticipated data distribution, the prompt-engineering process, and the trust placed in the foundation model. We derive the foundation prior as an exponential-tilted, generalized Bayesian update of the user's primitive prior, where a trust parameter governs the weight assigned to synthetic data. We then show how synthetic data and the associated foundation prior can be incorporated into standard statistical and econometric workflows, and discuss their use in applications such as refining complex models, informing latent constructs, guiding experimental design, and augmenting random-coefficient and partially linear specifications. By treating generative outputs as structured, explicitly subjective priors rather than as empirical observations, the framework offers a principled way to harness foundation models in empirical work while avoiding the conflation of synthetic ”facts” with real data.

We introduce Smart Bayes, a new classification framework that bridges generative and discriminative modeling by integrating likelihood-ratio-based generative features into a logistic-regression-style discriminative classifier. From the generative perspective, Smart Bayes relaxes the fixed unit weights of Naive Bayes by allowing data-driven coefficients on density-ratio features. From a discriminative perspective, it constructs transformed inputs as marginal log-density ratios that explicitly quantify how much more likely each feature value is under one class than another, thereby providing predictors with stronger class separation than the raw covariates. To support this framework, we develop a spline-based estimator for univariate log-density ratios that is flexible, robust, and computationally efficient. Through extensive simulations and real-data studies, Smart Bayes often outperforms both logistic regression and Naive Bayes. Our results highlight the potential of hybrid approaches that exploit generative structure to enhance discriminative performance.

Multipurpose batch processes become increasingly popular in manufacturing industries since they adapt to low-volume, high-value products and shifting demands. These processes often operate in a dynamic environment, which faces disturbances such as processing delays and demand changes. To minimise long-term cost and system nervousness (i.e., disruptive changes to schedules), schedulers must design rescheduling strategies to address such disturbances effectively. Existing methods often assume complete look-ahead information over the scheduling horizon. This assumption contrasts with realistic situations where schedulers can only access incomplete look-ahead information. Sticking with existing methods may lead to suboptimal long-term costs and high-level system nervousness. In this work we propose a Bayesian dynamic scheduling method. Our method relies on learning a Bayesian Network from the probability distribution of disturbances. Specifically, the Bayesian Network represents how likely each operation will be impacted by disturbances. During the online execution, when new disturbances become observed, this method updates the posterior distribution and therefore guides the rescheduling strategy. We compare our method with the existing periodic rescheduling strategy (which generates new schedules from scratch at fixed intervals) on four benchmark problems. Computational results show that our method achieves statistically better long-term costs and system nervousness. In the theoretical aspect, we prove that if disturbances are mutually independent, the impact-quantifying variables inherently satisfy the independence assumptions required by Bayesian Networks. As an implication, practitioners can extend the method to other scheduling problems (such as job shop scheduling and continuous processes), given that they define the problem-specific dependencies between operations.

Accurate and reliable forecasting models are critical for guiding public health responses and policy decisions during pandemics such as COVID-19. Retrospective evaluation of model performance is essential for improving epidemic forecasting capabilities. In this study, we used COVID-19 wastewater data from CDC's National Wastewater Surveillance System to generate sequential weekly retrospective forecasts for the United States from March 2022 through September 2024, both at the national level and for four major regions (Northeast, Midwest, South, and West). We produced 133 weekly forecasts using 11 models, including ARIMA, generalized additive models (GAM), simple linear regression (SLR), Prophet, and the n-sub-epidemic framework (top-ranked, weighted-ensemble, and unweighted-ensemble variants). Forecast performance was assessed using mean absolute error (MAE), mean squared error (MSE), weighted interval score (WIS), and 95% prediction interval coverage. The n-sub-epidemic unweighted ensembles outperformed all other models at 3-4-week horizons, particularly at the national level and in the Midwest and West. ARIMA and GAM performed best at 1-2-week horizons in most regions, whereas Prophet and SLR consistently underperformed across regions and horizons. These findings highlight the value of region-specific modeling strategies and demonstrate the utility of the n-sub-epidemic framework for real-time outbreak forecasting using wastewater surveillance data.

We study the geometry of the fixed-rank core covariance manifold arising from the Kronecker-core decomposition of covariance matrices. As shown in Hoff, McCormack, and Zhang (2023), every covariance matrix $\Sigma$ of $p_1\times p_2$ matrix-variate data uniquely decomposes into a separable component $K$ and a core component $C$. Such a decomposition also exists for rank-$r$ $\Sigma$ if $p_1/p_2+p_2/p_1<r$, with $C$ sharing the same rank. If this core $C$ exhibits a partial-isotropy structure, then a partial-isotropy rank-$r$ core is a non-trivial convex combination of a rank-$r$ core and $I_p$ for $p:=p_1p_2$, where the weight on $I_p$ measures the deviation of $\Sigma$ from separability. This motivates studying the geometry of the space of rank-$r$ cores, $\mathcal{C}_{p_1,p_2,r}^+$. We show that $\mathcal{C}_{p_1,p_2,r}^+$ is a smooth manifold, except for a measure-zero subset associated with canonical decomposability. When $r=p$, $\mathcal{C}_{p_1,p_2}^{++}:=\mathcal{C}_{p_1,p_2,p}^+$ is itself a smooth manifold. The geometric properties, including smoothness of the positive definite cone via separability and the Riemannian gradient and Hessian operator relevant to $\mathcal{C}_{p_1,p_2,r}^+$, are also derived. As an application, we propose a partial-isotropy core shrinkage estimator for matrix-variate data.

Monitoring the safety of medical products is a core concern of contemporary pharmacovigilance. To support drug safety assessment, Spontaneous Reporting Systems (SRS) collect reports of suspected adverse events of approved medical products offering a critical resource for identifying potential safety concerns that may not emerge during clinical trials. Modern nonparametric empirical Bayes methods are flexible statistical approaches that can accurately identify and estimate the strength of the association between an adverse event and a drug from SRS data. However, there is currently no comprehensive and easily accessible implementation of these methods. Here, we introduce the R package pvEBayes, which implements a suite of nonparametric empirical Bayes methods for pharmacovigilance, along with post-processing tools and graphical summaries for streamlining the application of these methods. Detailed examples are provided to demonstrate the application of the package through analyses of two real-world SRS datasets curated from the publicly available FDA FAERS database.

Clinical trials assessing neurological treatment are challenging due to the diversity of brain function, and the difficulty in quantifying it. Traditional treatment studies in epilepsy use seizure frequency as the primary outcome measure, which may overlooking meaningful improvements in patients' quality of life. This paper introduces the Clinical Instrument for Measuring Patient Anecdotes in Clinical Trials (Clinical IMPACT), a novel tool designed to capture qualitative non-seizure improvement across neurological domains. The Clinical IMPACT incorporates open-ended inquiries that allow participants or caregivers to identify and select anecdotal evidence of their most significant treatment benefits. A blinded panel of experts ranks these anecdotes, facilitating a rigorous statistical analysis using the Wilcoxon Rank-Sum Test to detect treatment efficacy. The approach is resistant to type 1 error, yet comprehensive in its ability to capture real-world effects on quality of life. The potential of the Clinical IMPACT tool to enhance sensitivity while also providing qualitative insights that can inform patients, healthcare providers, and regulatory bodies about treatment effects makes it important to consider in any neurological trial. We describe how it can be used in epilepsy, and advocate for its inclusion as a key secondary endpoint to provide a perspective on non-seizure outcomes, which have previously been challenging to measure, let alone to interpret, even when the clinical trial is positive.

We consider a statistical model of a n-mode quantum Gaussian state which is shift invariant and also gauge invariant. Such models can be considered analogs of classical Gaussian stationary time series, parametrized by their spectral density. Defining an appropriate quantum spectral density as the parameter, we establish that the quantum Gaussian time series model is asymptotically equivalent to a classical nonlinear regression model given as a collection of independent geometric random variables. The asymptotic equivalence is established in the sense of the quantum Le Cam distance between statistical models (experiments). The geometric regression model has a further classical approximation as a certain Gaussian white noise model with a transformed quantum spectral density as signal. In this sense, the result is a quantum analog of the asymptotic equivalence of classical spectral density estimation and Gaussian white noise, which is known for Gaussian stationary time series. In a forthcoming version of this preprint, we will also identify a quantum analog of the periodogram and provide optimal parametric and nonparametric estimates of the quantum spectral density.

Results in epidemiology and social science often require the removal of confounding effects from measurements of the pairwise correlation of variables in survey data. This is typically accomplished by some variant of linear regression (e.g., “logistic" or “Cox proportional"). But, knowing whether all possible confounders have been identified, or are even visible (not latent), is in general impossible. Here, we exhibit two examples that frame the issue. The first example proposes a highly unlikely hypothesis on drug use, draws data from a large, respected survey, and succeeds in “proving" the implausible hypothesis, despite regressing out more than 20 confounding variables. The second constructs a “metamodel" in which a single (by hypothesis unmeasurable) latent variable affects many mutually correlated confounders. From simulations, we derive formulas for the magnitude of spurious association that persists even as increasing numbers of confounders are regressed out. The intent of these examples is for them to serve as cautionary tales.

The wavelet spectra is a common starting point for estimating the Hurst exponent of a self-similar signal using wavelet-based techniques. The decay of the $\log_2$ average energy of the detail wavelet coefficients as a function of the level of signal decomposition can be used to construct estimators for this parameter. In this paper, we expand on previous work which introduced the “dual" wavelet spectra, where decomposition levels are instead treated as a function of energy values, and propose a relationship between its slope and the Hurst exponent by inverting the standard wavelet spectra, thereby creating a new estimator. The effectiveness of this estimator and its sensitivity to several settings are demonstrated through a simulation study. Finally, we show how the technique performs as a feature extraction method by applying it to the task of detecting the presence of breast cancer in mammogram images. Dual spectra wavelet features had a statistically significant effect on the log-odds of Cancer.

Principal Component Analysis (PCA) and K-means constitute fundamental techniques in multivariate analysis. Although they are frequently applied independently or sequentially to cluster observations, the relationship between them, especially when K-means is used to cluster variables rather than observations, has been scarcely explored. This study seeks to address this gap by proposing an innovative method that analyzes the relationship between clusters of variables obtained by applying K-means on transposed data and the principal components of PCA. Our approach involves applying PCA to the original data and K-means to the transposed data set, where the original variables are converted into observations. The contribution of each variable cluster to each principal component is then quantified using measures based on variable loadings. This process provides a tool to explore and understand the clustering of variables and how such clusters contribute to the principal dimensions of variation identified by PCA.

Echocardiogram datasets enable training deep learning models to automate interpretation of cardiac ultrasound, thereby expanding access to accurate readings of diagnostically-useful images. However, the gender, sex, race, and ethnicity of the patients in these datasets are underreported and subgroup-specific predictive performance is unevaluated. These reporting deficiencies raise concerns about subgroup validity that must be studied and addressed before model deployment. In this paper, we show that current open echocardiogram datasets are unable to assuage subgroup validity concerns. We improve sociodemographic reporting for two datasets: TMED-2 and MIMIC-IV-ECHO. Analysis of six open datasets reveals no consideration of gender-diverse patients and insufficient patient counts for many racial and ethnic groups. We further perform an exploratory subgroup analysis of two published aortic stenosis detection models on TMED-2. We find insufficient evidence for subgroup validity for sex, racial, and ethnic subgroups. Our findings highlight that more data for underrepresented subgroups, improved demographic reporting, and subgroup-focused analyses are needed to prove subgroup validity in future work.

We study the multi-objective linear contextual bandit problem, where multiple possible conflicting objectives must be optimized simultaneously. We propose MOL-TS, the first Thompson Sampling algorithm with Pareto regret guarantees for this problem. Unlike standard approaches that compute an empirical Pareto front each round, MOL-TS samples parameters across objectives and efficiently selects an arm from a novel effective Pareto front, which accounts for repeated selections over time. Our analysis shows that MOL-TS achieves a worst-case Pareto regret bound of $\widetilde{O}(d^{3/2}\sqrt{T})$, where $d$ is the dimension of the feature vectors, $T$ is the total number of rounds, matching the best known order for randomized linear bandit algorithms for single objective. Empirical results confirm the benefits of our proposed approach, demonstrating improved regret minimization and strong multi-objective performance.

We address the problem of causal effect estimation in the presence of hidden confounders using nonparametric instrumental variable (IV) regression. An established approach is to use estimators based on learned spectral features, that is, features spanning the top singular subspaces of the operator linking treatments to instruments. While powerful, such features are agnostic to the outcome variable. Consequently, the method can fail when the true causal function is poorly represented by these dominant singular functions. To mitigate, we introduce Augmented Spectral Feature Learning, a framework that makes the feature learning process outcome-aware. Our method learns features by minimizing a novel contrastive loss derived from an augmented operator that incorporates information from the outcome. By learning these task-specific features, our approach remains effective even under spectral misalignment. We provide a theoretical analysis of this framework and validate our approach on challenging benchmarks.

Evaluating rare-event forecasts is challenging because standard metrics collapse as event prevalence declines. Measures such as F1-score, AUPRC, MCC, and accuracy induce degenerate thresholds – converging to zero or one – and their values become dominated by class imbalance rather than tail discrimination. We develop a family of rare-event-stable (RES) metrics whose optimal thresholds remain strictly interior as the event probability approaches zero, ensuring coherent decision rules under extreme rarity. Simulations spanning event probabilities from 0.01 down to one in a million show that RES metrics maintain stable thresholds, consistent model rankings, and near-complete prevalence invariance, whereas traditional metrics exhibit statistically significant threshold drift and structural collapse. A credit-default application confirms these results: RES metrics yield interpretable probability-of-default cutoffs (4-9%) and remain robust under subsampling, while classical metrics fail operationally. The RES framework provides a principled, prevalence-invariant basis for evaluating extreme-risk forecasts.

Grouped multiple testing arises in many high-dimensional applications, where hypotheses may differ across groups and valid p-values are unavailable. Applying competition-based procedures separately within groups can inflate the global false discovery rate (FDR), while conservative corrections may substantially reduce power. We extend the competition test framework to the grouped setting and propose two procedures. The group competition (GC) test provides finite-sample FDR bounds under arbitrary inter-group dependence. For homogeneous groups, the required correction has logarithmic order in the number of groups, and this order is unavoidable for the GC test. The consistent group competition (CGC) filter instead couples group-specific thresholds through a single global normalization. Under the stated basic-independence conditions, it controls the FDR in finite samples without a correction that grows with the number of groups; under arbitrary inter-group dependence, the same guarantee is retained with a per-group correction that depends on the symmetry parameter of each group. We further develop a parameterized-threshold analysis to establish asymptotic guarantees and to characterize when the per-group correction can be dropped, so that a simple filter still controls the FDR approximately. Simulation studies show that the GC test and the CGC filter can improve power over pooled competition testing in heterogeneous settings, whereas the pooled test may remain more powerful under homogeneity. A peptide-spectrum matching benchmark involving post-translational modifications illustrates the practical relevance of the proposed methods.

Large spatial datasets with non-Gaussian responses are increasingly common in environmental monitoring, ecology, and remote sensing, yet scalable Bayesian inference for such data remains challenging. Markov chain Monte Carlo methods are often prohibitive for large datasets, and existing variational Bayes methods rely on conjugacy or strong approximations that limit their applicability and can underestimate posterior variances. A scalable variational framework that incorporates semi-implicit variational inference (SIVI) with basis representations of spatial generalized linear mixed models, which may not have conjugacy, is proposed. The proposed framework accommodates gamma, negative binomial, Poisson, Bernoulli, and Gaussian responses on continuous spatial domains. Across 20 simulation scenarios with 50,000 locations, SIVI achieves predictive accuracy and posterior distributions comparable to Metropolis–Hastings and Hamiltonian Monte Carlo while providing notable computational speedups. Applications to remotely-sensed land surface temperature and blue jay abundance further demonstrate the utility of the approach for large non-Gaussian spatial datasets.

This paper introduces ARCADIA, an agentic AI framework for causal discovery that integrates large-language-model reasoning with statistical diagnostics to construct valid, temporally coherent causal structures. Unlike traditional algorithms, ARCADIA iteratively refines candidate DAGs through constraint-guided prompting and causal-validity feedback, leading to stable and interpretable models for real-world high-stakes domains. Experiments on corporate bankruptcy data show that ARCADIA produces more reliable causal graphs than NOTEARS, GOLEM, and DirectLiNGAM while offering a fully explainable, intervention-ready pipeline. The framework advances AI by demonstrating how agentic LLMs can participate in autonomous scientific modeling and structured causal inference.

Differential privacy is increasingly formalized through the lens of hypothesis testing via the robust and interpretable $f$-DP framework, where privacy guarantees are encoded by a baseline Blackwell trade-off function $f_{\infty} = T(P_{\infty}, Q_{\infty})$ involving a pair of distributions $(P_{\infty}, Q_{\infty})$. The problem of choosing the right privacy metric in practice leads to a central question: what is a statistically appropriate baseline $f_{\infty}$ given some prior modeling assumptions? The special case of Gaussian differential privacy (GDP) showed that, under compositions of nearly perfect mechanisms, these trade-off functions exhibit a central limit behavior with a Gaussian limit experiment. Inspired by Le Cam's theory of limits of statistical experiments, we answer this question in full generality in an infinitely divisible setting. We show that suitable composition experiments $(P_n^{\otimes n}, Q_n^{\otimes n})$ converge to a binary limit experiment $(P_{\infty}, Q_{\infty})$ whose log-likelihood ratio $L = \log(dQ_{\infty} / dP_{\infty})$ is infinitely divisible under $P_{\infty}$. Thus any limiting trade-off function $f_{\infty}$ is determined by an infinitely divisible law $P_{\infty}$, characterized by its Levy–Khintchine triplet, and its Esscher tilt defined by $dQ_{\infty}(x) = e^{x} dP_{\infty}(x)$. This characterizes all limiting baseline trade-off functions $f_{\infty}$ arising from compositions of nearly perfect differentially private mechanisms. Our framework recovers GDP as the purely Gaussian case and yields explicit non-Gaussian limits, including Poisson examples. It also positively resolves the empirical $s^2 = 2k$ phenomenon observed in the GDP paper and provides an optimal mechanism for count statistics achieving asymmetric Poisson differential privacy.