Papers

The Donsker–Varadhan formula characterizes the ordinary Bayesian posterior as the solution of an unrestricted $\mathsf{KL}$-regularized variational problem. Generalized variational inference replaces this regularizer by other divergences, but the resulting measure-valued optimization problem is often studied only after restriction to a parametric variational family. This paper studies the unrestricted measure-level problem. Given a measurable space $(\mathcal{Z},\mathfrak{Z})$, a prior probability measure $P$, a measurable loss $\ell:\mathcal{Z}\to(-\infty,\infty]$, a regularization strength $\alpha>0$, and a divergence $\mathsf{D}(Q\Vert P)$, we seek probability measures in \[ \underset{Q\in\mathcal{P}(\mathcal{Z})}{\mathrm{arg\,min}}\left\{\int_{\mathcal{Z}} \ell\,\mathrm{d}Q+\alpha\mathsf{D}(Q\Vert P)\right\}. \] For $f$-divergence penalties we derive a scalar inverse-gradient density formula and a one-dimensional dual identity; the Kullback–Leibler, Cressie–Read, and squared-Hellinger problems are treated as examples. Reverse $f$-divergences and mixed forward/reverse Kullback–Leibler penalties follow from the same separable integral principle. For Bregman divergences between densities we obtain a density-space solution with a scalar mass multiplier, including least-squares, density-power, and Burg/Itakura–Saito examples. For Rényi penalties of order $r>1$ we derive a normalized truncated-power characterization and a threshold equation for every global optimizer. Finite model-weight formulas and simple conjugate Bayesian model illustrations show how these closed forms are realized in practice and differ from the traditional solutions.

Community detection is a fundamental problem in network analysis. While many existing methods focus on homogeneous networks, real world networks are often heterogeneous, involving multiple node types and interaction mechanisms. In addition, node specific covariates frequently provide valuable information about the underlying community structure. Existing methodologies typically account for either network heterogeneity or covariate information, but seldom both simultaneously. In this paper, we propose a covariate assisted spectral clustering framework for heterogeneous networks that jointly utilizes network connectivity in a heterogeneous setting and node level covariates. The proposed method extends covariate assisted spectral clustering to heterogeneous settings and operates directly on the heterogeneous network without relying on projection based simplifications. Under a heterogeneous node contextualized stochastic blockmodel, we establish theoretical guarantees for the proposed procedure, including concentration results, eigenspace perturbation bounds, and an explicit upper bound on the misclustering rate. Simulation studies demonstrate that incorporating covariate information substantially improves community recovery and consistently outperforms several benchmark methods. We further apply the proposed framework to United Nations General Assembly voting data, where it reveals meaningful geopolitical structures by combining voting interactions with auxiliary covariate information.

Existing approaches to e-scooter mobility hub planning lack city-type-specific causal evidence. Demand models are typically correlational, built on proprietary trip data, and do not distinguish how driver profiles vary across urban typologies. This paper presents a three-phase agentic AI framework that constructs a Causal Template Library from public GBFS data across 29 German cities, encoding which environmental features causally drive hotspot demand for each combination of city type (large, university, industrial, hilly) and cluster type (core, peripheral). A large language model (LLM) orchestrated causal discovery pipeline adapts algorithm selection to local data conditions across 57 city-cluster units. The library reveals systematic variation. Core demand is driven by activity access and transit proximity, while peripheral demand responds to built form, with city-type-specific patterns supporting transferable siting templates. A planning tool built on the library scores candidate sites, calibrates infrastructure recommendations to local demographics, and generates practitioner-ready reports. In Heilbronn, Germany, two hub sites informed by the framework's causal evidence are currently under construction, illustrating how the outputs can support real-world siting decisions.

We provide an upper-bound for the Lévy concentration function: $$ Q_{S}(\varepsilon):= \sup_{x \in\mathbb{R}}\mathbb{P} (x < S \leq x+\varepsilon) $$ where $S$ is a weighted sum of noncentral chi-square random variables: $$ S:= \sum_{k=1}^\infty \lambda_k (Z_k^2 - 1) + \mu_kZ_k $$ Here, $\{Z_k\}_{k=1}^\infty$ is a sequence of independent standard Gaussian random variables and $\{\lambda_k\}_{k=1}^\infty, \{\mu_k\}_{k=1}^\infty$ are real valued, square summable sequences. Random variables of this type often appear as limiting distributions of second order U-statistics. Our bound is adaptive, in that it recovers (up to constant factors) Gaussian type concentration function estimates if $\|\lambda\|_2$ is negligible compared to $\|\mu\|_2$ and chi-square estimates if $\|\mu\|_{2}$ is negligible compared to $\|\lambda\|_2$. Our bound generalizes existing bounds in various ways. In particular, we make no assumptions regarding the number of nonzero $|\lambda_k|$ or the size of the minimal $|\lambda_k|$, nor do we make any assumptions on the signs of $\lambda_k$. Finally, we apply our bound to some examples of interest, specifically quadratic forms that arise in limit theorems for second-order U-statistics.

We develop a general probabilistic framework for analyzing propagation of chaos in transport ensemble filters (TEFs), a broad class of interacting particle systems that are used to approximate the sequence of state distributions in hidden Markov models given a history of observations. This class of transport-based filtering algorithms includes the widely used ensemble Kalman filter (EnKF), based on affine updates at each filtering step, as well as the ensemble stochastic map filter (EnSMF), which employs nonlinear updates. For this class, we identify the limiting mean-field dynamics. We then establish non-asymptotic, high-probability, pathwise Wasserstein convergence of the interacting particle system to an i.i.d. ensemble drawn from this mean-field limit at the Monte Carlo rate. Convergence to the mean-field law itself follows with the usual dimension-dependent empirical Wasserstein rate. The proof combines a synchronous coupling construction with stability of moments and tails under conditioning, together with quantitative estimates for the propagation of the underlying dynamics through the interacting particle system. Applying our theory to both the EnKF and the EnSMF yields the first non-asymptotic, high-probability convergence guarantees for TEFs.

As black-box models become foundational to modern research, ensuring their stability is paramount for the realization of trustworthy artificial intelligence. The inherent diversity of inputs - ranging from structured Gaussian distributions to complex data with unknown structures - poses a significant challenge: how to stabilize black-box outputs while effectively leveraging available prior information. This paper introduces a task-oriented randomization methodology that adaptively tailors its strategy to the underlying generative mechanisms of the input data, specifically addressing unstructured complexities. A comprehensive suite of stability guarantees is proposed. Beyond establishing rigorous theoretical foundations for stability, the research provides a detailed analysis of the intrinsic trade-off between stability and exploration. Motivated by the architecture of Large Language Models, the framework is further extended to top-k ranking problems. The validity and effectiveness of the proposal are demonstrated through extensive numerical simulations and applications to the real-world dataset.

Textual data are often collected alongside structured response variables, but prediction and interpretation are commonly treated as separate tasks. This paper studies rating prediction as an initial case of interpretable text-response modelling, where the aim is to learn textual representations that are both semantically meaningful and aligned with an external response. We propose a joint non-negative matrix factorisation and binomial regression model, in which the document-topic representation is learned from both text reconstruction and rating prediction. Simulation experiments and a real-world review dataset show that the model can recover stable response-relevant textual signals and achieve competitive performance against linear and ridge regression baselines. The framework provides a practical step towards interpretable modelling of text-linked outcomes, with potential extensions to other response types beyond bounded ratings.

Social disadvantage is associated with multimorbidity, but the pathways linking social conditions to disease burden remain poorly understood. We developed an AI-driven multimodal mediation framework that integrates socioeconomic, psychosocial, clinical, laboratory, behavioral, and genomic data from the All of Us Research Program. Modality-specific variational autoencoders were used to derive latent representations of each data domain, and mediation analyses were subsequently performed in latent space to evaluate indirect associations between socioeconomic disadvantage, psychosocial factors, and multimorbidity. The final analytic cohort included 20,804 participants with complete multimodal data. Across 800 exposure–mediator–outcome combinations, mediation signals were concentrated within a small number of latent dimensions. The strongest indirect association linked a socioeconomic disadvantage dimension, a psychosocial vulnerability dimension, and a cardiometabolic multimorbidity dimension (NIE = 0.002517). The psychosocial dimension was characterized by poorer mental health, greater loneliness, lower social well-being, and lower health literacy, whereas the outcome dimension was associated with hypertension, diabetes, hyperlipidemia, obesity, chronic kidney disease, and heart disease. Bootstrap analyses supported the stability of the leading pathway. These findings suggest that psychosocial vulnerability was strongly represented in the dominant latent pathway linking socioeconomic disadvantage and cardiometabolic multimorbidity. More broadly, the proposed framework illustrates how AI-based representation learning can be used to investigate complex relationships across high-dimensional multimodal health data.

Tabular prediction is central to a wide range of real-world applications. Tabular data typically contain heterogeneous features as well as rich and complex relational information that can imply a latent structural manifold. Hyperbolic geometry can help capture complex structural relations in data. However, most existing tabular prediction models are constructed and optimized in Euclidean space. How to incorporate hyperbolic geometry into supervised tabular learning remains underexplored. We propose HTNN, a supervised hyperbolic manifold constrained tabular neural network for tabular prediction. HTNN consists of a hyperbolic feature-value representation layer for heterogeneous categorical and numerical features, followed by a conventional MLP predictor. HTNN employs a geometry-aware training and geometry-free inference optimization framework. The geometry-aware training allows hyperbolic geometry to shape the latent representation learning of heterogeneous feature values. After training, the latent hyperbolic representations can be converted into ordinary Euclidean space for efficient geometry-free inference. We conducted extensive experiments on the TALENT benchmark. HTNN ranks first among 36 methods on 200 classification datasets and third among 34 methods on 100 regression datasets. Experimental results show that the proposed hyperbolic manifold constrained tabular neural network is effective.

Learning effective policies for adaptive data acquisition remains challenging: posterior-based methods rely on surrogate models and posterior approximations that can be misspecified or biased, while direct policy-learning methods map from historical observations and fail to exploit available model representations, making learning harder. We introduce policy learning with belief representations (POLAR), based on the insight that optimal data acquisition depends on the observation history only through a sufficient belief state. Specifically, POLAR decouples representation learning from policy learning by leveraging pretrained predictive foundation models as belief-state encoders, training a policy head on top of their representations. This yields a simple, unified amortised policy learning framework for Bayesian experimental design, Bayesian optimisation, and active learning, differing only in the task-specific utility used to train the policy. Empirically, we find that POLAR outperforms state-of-the-art amortised methods across diverse tasks while requiring far fewer training samples, demonstrating a significant step in the scalability and efficiency of amortised data acquisition.

We introduce a PAC framework in which the learner can access sampling oracles both before and at decision time. Sample complexity is measured by a pair $(n,m)$, where $n$ is the learning budget spent before a query is known and $m$ is the additional sampling budget per query. We demonstrate its relevance in discounted Markov decision processes with exogenous i.i.d.\ contexts revealed before acting. Contexts may affect both rewards and transitions but remain uncontrolled by the agent. The learner can sample the unknown context distribution and the transition kernel. We study policy evaluation (PE), best-value estimation (BVE), and best-policy extraction (BPE). When rewards and transitions are known, a variance-reduced algorithm solves all three tasks with sample complexity $\bigl(\widetilde O((1-\gamma)^{-3}\varepsilon^{-2}),0\bigr)$, which is minimax optimal up to logarithmic factors. Let $\mathcal{X}$ be the controlled state space. When transitions are also unknown, we give a PE algorithm with complexity $\bigl(\widetilde O(|\mathcal X|(1-\gamma)^{-3}\varepsilon^{-2}), \widetilde O((1-\gamma)^{-2}\varepsilon^{-2})\bigr)$ and matching lower bounds at this budget pair. For BVE and BPE, we give an algorithm with a common offline budget $\widetilde O(|\mathcal X|^2|\mathcal A|(1-\gamma)^{-4}\varepsilon^{-2})$ and respective query costs $\widetilde O(|\mathcal A|(1-\gamma)^{-2}\varepsilon^{-2})$ and $\widetilde O(|\mathcal A|(1-\gamma)^{-3}\varepsilon^{-2})$. Importantly, all bounds are independent of the context-space cardinality.

Sampling from an unnormalized target density by reversing an Ornstein-Uhlenbeck diffusion requires the score of each noise-perturbed marginal law. Two exact identities are available: Tweedie's identity and a target-score identity, each yielding unbiased finite-reference score estimators for the OU-marginal score. Score estimators induced by scalar blends of Tweedie and TSI score estimators can reduce variance, but they are too rigid for singular or strongly anisotropic targets. We formulate blended score estimation as a conditional risk-minimization problem over matrix valued blending coefficients, referred to as gates. Our central result is to show the optimal matrix valued gate for blended score estimation is given \[ G_\star(y,t) = \alpha_t^2 \left(\alpha_t^2 I_d + \gamma_t\, \mathbb{E}[H_0(X_0)\mid Y_t=y] \right)^{-1}, \qquad H_0=-\nabla^2\log p_0 .\] Here $\alpha_t = e^{-t}$ and $\gamma_t = 1-e^{-2t}$ are the OU coefficients, and the conditional expectation is under the OU posterior of $X_0$ given $Y_t=y$. We call this formula the Laplace-Fisher Gate Identity (\LFGI{}). Because the Tweedie-TSI disagreement has conditional mean zero, the gate changes the score-estimator variance but not its expected value. We derive the variance-optimal matrix gate, record the Gaussian special case, and establish finite-reference consistency and stability bounds for estimating the gate from weighted reference samples. We then use the finite-reference LFGI score estimator for normalized density evaluation in Bayesian inverse problems. In regimes where MCMC pilot samples and derivative information are already available, LFGI uses those byproducts to construct a normalized surrogate for the posterior density. The resulting surrogate supplies information that the MCMC samples alone do not provide: posterior-energy evaluation, model-evidence estimation, and downstream density-based diagnostics. On a PDE-constrained inverse-problem benchmark, the LFGI surrogate improves posterior-density calibration and sampling diagnostics relative to the other tested score-estimator classes. Experiments using LFGI with known model evidence check absolute evidence calibration in both Gaussian and non-Gaussian settings.

Mobile and wearable technologies enable the delivery of just-in-time adaptive interventions (JITAIs) – interventions that adapt treatment delivery to an individual's rapidly changing internal state and context in real-time, real-world settings. Estimating optimal JITAIs, however, remains challenging because these studies often involve dozens of decision points per individual, and existing methods can produce unstable and irregular estimators with substantial bias and slow convergence rates. Advanced reinforcement learning approaches may be difficult to interpret and often target proximal, discounted outcomes rather than the distal end-of-study outcomes that define long-term success in many behavioral and clinical studies. To address these challenges, we develop a nonparametrically efficient estimator of the regimen-response curve for distal outcomes under a class of stochastic policies and introduce a data-adaptive tilting procedure to stabilize estimation in settings with many decision points. We show that the estimated regimen-response curve converges weakly to a Gaussian process, enabling simultaneous confidence bands, and we derive asymptotic theory for the optimizer of the curve, thereby enabling inference for the learned optimal stochastic policy. These developments provide a unified framework for estimation, inference, and optimization of stochastic JITAIs for distal outcomes.

Many modern Language Model (LM) pipelines return an averaged model, such as an exponential moving average of the training iterates, rather than the final iterate itself. This raises a fundamental question: given that we will return an iterate average, how should we change training to improve the performance of this average? We study this question by formulating optimizer design for the iterate-average estimator as an optimal-control problem. In a continuous-time stochastic quadratic model, we solve for the control strategy that minimizes the error of the returned average subject to a penalty on the size of the intervention. A practical approximation to this controller yields PACE, a lightweight wrapper around AdamW that pulls the live weights toward their exponential moving average with a clipped, per-coordinate control strength. We prove that a stylized version of PACE converges at the standard stochastic convex optimization rate, up to a factor depending on the averaging rule, while in the quadratic setting it can strictly improve the limiting squared error of the iterate-average estimator and can do so by an arbitrarily large factor on some instances. Empirically, our results suggest that PACE improves over AdamW and EMA-evaluated AdamW in supervised fine-tuning of 1-2B parameter LMs and in GPT-2 pretraining on FineWeb for a wide range of learning rates, decay schedules, and other hyperparameters.

Small area estimation and disease mapping increasingly rely on areal data where reporting boundaries change over time. We develop a computationally efficient spatio-temporal disaggregation method to recover high-resolution risk surfaces from observed counts under changing boundaries. Our approach extends the spatially aggregated log-Gaussian Cox process and uses the Extended Latent Gaussian Model framework for fast approximate posterior inference. We replace standard lognormal polygon-specific effects with gamma-distributed overdispersion which yields a marginal negative binomial likelihood, and removes one latent variable per polygon-time pair. We illustrate the approach by mapping mortality risk across shifting NUTS-3 boundaries in Belgium and the Netherlands. For the purpose of dissemination we use Codex to leverage the methodology presented in this paper for the analysis of a separate data set concerning the city of Manchester. The methodology is implemented in the open-source R package DAST.

Rank aggregation combines information from ordered lists ranking items by preference. Classical parametric models for such data, including the Mallows and Plackett-Luce models, assume the orders concentrate around one or more complete consensus rankings. Recent work relaxes the total-order assumption by allowing the consensus structure to be a partial order (poset), allowing for incomparabilities in preferences. However, in many applications preference data exhibit group structure. We introduce hierarchical partial order (HPO) models, which extend poset-based models to accommodate grouped data through a hierarchy of latent posets. This framework, which parallels mixture model extensions of the Mallows and Plackett-Luce models, enables principled sharing of information across groups while preserving partial-order structure. We show that the Plackett-Luce model and its hierarchical variants are special cases of HPO-models. We develop a hierarchical clustering extension (HCPO) for unsupervised clustering in settings where group labels are unknown. Bayesian inference for the latent poset hierarchy is performed using Markov chain Monte Carlo methods. Experiments on synthetic and real-world datasets, including pairwise acoustic preference data and LLM agent traces, demonstrate that the proposed HPO and HCPO models outperform existing approaches in both predictive performance and structural interpretability.

We prove a single algebraic mixed coincidence identity that unifies a broad swath of information-theoretic variational results. For any family of priors $\{\pi_i\}$ and real exponents $\{ \alpha_i \}$, the log of the mixed count $E_{x\sim\nu}\!\left[\prod_{i=1}^W \pi_i^{\alpha_i}(x)\right]$ is simultaneously a Boltzmann coincidence weight, an exponential-family normalizer, a maximum-entropy value, and a KL-barycenter optimum. The identity yields a unified derivation of classical cornerstones of information theory: concentration of empirical distributions (Sanov-type decompositions and Gibbs conditioning), hypothesis-testing error exponents (Chernoff information and its multi-way analogue), change-of-measure inequalities (Donsker-Varadhan and PAC-Bayes), and laws governing rare-pattern coincidences (Erdos-Renyi run-length, iterative guesswork, rate-distortion, and birthday thresholds). Each is recovered as a specialization of the same algebraic equality. It strictly generalizes the classical Renyi entropy and divergence variational formulas (one and two priors respectively) to a $W$-prior simplex, and holds for unnormalized and continuum-indexed priors. Among its consequences are an exact multi-prior PAC-Bayes penalty that subtracts an explicit "coincidence bonus" from the usual single-prior posterior penalty, and the asymptotic MAP error exponent for $W$-ary hypothesis testing as an edge-restricted simplex optimum. We demonstrate the calculus at scale on two large alphabets encoding richly modeled sequential languages: on language-model next-token predictives where we recover contrastive decoding, and on human genomic regulatory sequence where it separates correlated from diverse prior families along a sliding-window trace.

Optimal transport (OT) has become a central language for comparing probability measures, but exact balanced OT is often both too rigid for data with missing, created, or destroyed mass and subject to unfavorable high-dimensional sample complexity. Entropic regularization and unbalanced relaxations address these limitations in complementary ways. Entropy smooths the geometry, improves statistical behavior, and enables fast Sinkhorn-type algorithms, while unbalanced marginal penalties replace hard conservation constraints by divergence terms adapted to noisy empirical data. This paper studies the sample complexity of entropic unbalanced OT at the level of the optimal coupling, rather than only the scalar transport value. We develop a translation-invariant dual formulation, prove compactness and strong convexity properties for the intrinsic dual variables, and convert these geometric estimates into high-probability finite-sample bounds for empirical couplings. The results clarify why regularization is a practical necessity in machine learning applications: it softens the curse of dimensionality, reduces the number of samples needed for stable transport estimation, and keeps the resulting estimators compatible with scalable Sinkhorn-type solvers.

Modeling and sampling from the underlying distribution of asynchronous event sequences are crucial in various real-world applications, including social networks, medical diagnosis, and financial transactions. Existing autoregressive methods suffer from error accumulation during multi-step generation, while non-autoregressive diffusion methods are typically limited to fixed-length output sequences. In this paper, we propose Latent Block-Diffusion Temporal Point Processes (LBDTPP), a novel semi-autoregressive TPP framework that introduces a latent block diffusion mechanism for high-quality and variable-length event sequence generation. The core idea is to define an autoregressive probability distribution over event blocks in latent space and perform Gaussian diffusion within each block. By sequentially generating blocks while simultaneously sampling events in each block, LBDTPP preserves the length flexibility of autoregressive TPPs and inherits the parallel high-quality generation capability of diffusion models. Theoretically, we derive Wasserstein error bounds showing that, under suitable local approximation and prefix-stability assumptions, block-wise generation can reduce error accumulation compared with event-wise autoregressive generation. Extensive experiments on six real-world benchmark datasets demonstrate that LBDTPP outperforms state-of-the-art TPP baselines in both unconditional and conditional generation tasks. Further empirical analyses verify the benefits of latent-space diffusion and block-wise generation, and reveal the trade-off between generation quality and block size. Our code is available at https://github.com/Zh-Shuai/LBDTPP.

We study linear TD(0) under Markovian sampling, where data are generated along a single trajectory. We provide high-probability guarantees for a plain unprojected TD(0) algorithm with Polyak-Ruppert (PR) averaging, using a single stepsize schedule $\eta_t \propto \frac{1}{\tau_{\mathrm{mix}}\log(t)\sqrt{t}}$ that depends on the mixing time but requires no prior knowledge of the curvature parameter $\omega$. Our first result shows that such a choice of the stepsize guarantees that the TD(0) iterates are automatically and uniformly bounded with high probability, without projections and without any stability argument based on $\omega$. Building on this result, we establish a simultaneous high-probability convergence guarantee for the PR average: the same stepsize yields both a robust curvature-free $\widetilde{\mathcal{O}}\!\left(\frac{\tau_{\mathrm{mix}}}{\sqrt{T}}\right)$ rate and a fast curvature-dependent $\widetilde{\mathcal{O}}\!\left(\frac{\tau_{\mathrm{mix}}^2}{\omega T}\right)$rate, with the bound taking the minimum of the two. The core technical ingredient is a Poisson-equation toolkit for geometrically mixing Markov chains, which decomposes Markov noise into a martingale term plus a controlled remainder and enables a new self-bounding inductive argument for pathwise stability.

Modern online experimentation platforms produce data continuously, yet practitioners often apply Fixed Horizon Testing (FHT) while repeatedly peeking at results, inflating Type I error. Popular always-valid methods permit early efficacy stopping but do not natively combine futility stopping, business-relevant launch thresholds, and Type II error control. We adapt Wald's Sequential Probability Ratio Test (SPRT) to meet these needs. SPRT-z uses a large-sample normal approximation to Hajnal's sequential t-test, reducing computation and enabling two Brownian methods. Scale-Free Horizon Calibration (SFHC) sets a maximum sample size to approximately recover fixed-horizon power under discrete monitoring and futility stopping. A Brownian median unbiased estimator (MUE) and confidence intervals are defined for every stopping outcome. A validity analysis provides exact Gaussian-oracle accounting for the deployed stopping rule and a qualified fixed-look approximation for stationary finite-variance outcomes. In simulations with identically distributed arms, SPRT-z kept Type I error at or below nominal, tracked fixed-horizon power within a few percentage points except under the coarsest monitoring, and reduced expected sample size in all but one underpowered setting. The MUE reduced early-stopping bias, with close-to-nominal interval coverage in most settings; we report exceptions and operating rules. Stress tests identify outcome shapes and monitoring schedules that require separate error and power calibration. Across production experiments at Atlassian, median sample saving relative to fixed-horizon plans was 45.5%. More than 75% of decisions came from the futility boundary, supporting its practical importance.

This paper introduces a new econometric framework for modeling social interactions with heterogeneous peer responses, addressing endogenous link formation. Our Selection-corrected Heterogeneous Spatial Autoregressive (SCHSAR) approach jointly models link formation and outcome determination. We incorporate a finite mixture structure to capture heterogeneity in peer effects and account for unobserved individual-specific factors driving both network formation and outcome equations, addressing network endogeneity for credible estimation of heterogeneous spillover effects. We propose a fully Bayesian data augmentation approach for estimation and inference, overcoming challenges posed to standard likelihood-based methods. A simulation study validates our approach. Our empirical application to an innovation network among U.S. firms reveals significant positive, yet heterogeneous, peer effects on corporate R&D investments, after accounting for endogenous network formation. The findings highlight varying firm behaviors in response to exogenous R&D policy shocks and and quantify firm-level direct and spillover effects, offering valuable insights for evidence-based and targeted policy design.

Agentic data analysis systems produce rich outputs, including code, numerical results, and verbal diagnostics. This makes them more challenging to evaluate than single-turn LLM responses. It is therefore necessary to distinguish genuine disagreement between an agent's output and a ground-truth answer from grading artifacts. We investigate how reliably automated graders assess such a system and what strategies improve grading quality by applying LAMBDA, a multi-agent data-analysis system, on 153 numerical QRData tasks from DSGym. We develop and evaluate a three-layer human-AI grading cascade: strict regex matching, LLM-based lenient grading, and snippet-based human inspection, which combines non-GenAI and GenAI strategies with different failure profiles. Both automated graders achieve 100% observed precision (0/70 false positives). The lenient grader's recall is 97% against human labels. A keyword-anchored extraction pipeline raises the strict grader's recall by 60 percentage points over a last-number heuristic; the lenient grader is architecturally parser-independent. An iterative nudge mechanism raises grading run success from 36% to 97% and lenient-pass rates from 16% to 46%; comparing nudging with and without original-question re-injection shows that re-injection offers no benefit, confirming the nudge as an answer template cue. We further observe in this case study that variable type is the task metadata field most consistently associated with grading pipeline dynamics and observed outcome grades.

This paper introduces a Riemannian autoregressive (R-AR) model of order one for manifold-valued time series. The model is specified through an autoregressive process in the tangent space at a reference point, which is mapped to the manifold via the exponential map. It is characterised by two parameters: a reference point $\mu$, representing the intrinsic central tendency through the Fréchet mean, and an autoregressive parameter $\phi$, governing the dependence structure and stationarity properties. When these parameters are unknown, their estimation introduces geometric and probabilistic challenges arising from the intrinsic estimation of $\mu$ and from estimating $\phi$ as the base point varies over the manifold and the associated tangent spaces change. The estimation theory is developed using a strong law of large numbers for the sample Fréchet mean set of ergodic Markov chains, together with ergodic arguments for the strong consistency of the autoregressive parameter estimator. The framework is validated through numerical simulations in the hyperbolic plane and an application to aerosol size distributions on the Fisher-Rao manifold, demonstrating how the proposed model can characterise mean-reverting dynamics in nonlinear geometries.

This paper presents general strong duality results when testing hypotheses by betting against them. A bet is an e-variable for a composite null hypothesis $\Pcal$: a nonnegative random variable $X$ whose expected value is at most one under every $P \in \mathcal P$. Following Kelly, Breiman, Cover, Shafer, and Grunwald et al. (2024), we study a natural minimax log-optimality criterion: given a composite alternative $\Qcal$, we characterize the “GROW value” $\sup_{X} \inf_{Q} \E_{Q}[\log X]$. This paper generalizes the results of Larsson et al. (2025) from (arbitrary $\mathcal P$ and) simple $\mathcal Q$ to arbitrary $\mathcal Q$. We prove that there always exists a minimizing information-projection pair between the weak-$*$ closures of the convex hulls of arbitrary $\mathcal P$ and $\mathcal Q$, and show that the GROW value for bounded e-variables always equals their relative entropy. We also prove a similarly general strong duality for the REGROW criterion with bounded e-variables and arbitrary bounded offsets. Under various assumptions our results extend to unbounded e-variables, and examples show that without any assumptions such extensions fail. Our results are analogous to those in Larsson et al. (2026), swapping tests for bounded e-variables, minimax risk for the GROW criterion, and total variation for relative entropy.

Concentration inequalities for sample covariance matrices are fundamental tools in high-dimensional probability. Classical results typically assume that the selected random vectors are independent of the selection rule. In this paper, we study spectral concentration for sample covariance matrices formed from arbitrary, possibly data-dependent subsets of i.i.d. random vectors. Such data-dependent selection destroys the usual independence structure and makes standard covariance concentration bounds inapplicable. For i.i.d. Gaussian random vectors, we prove high-probability lower and upper bounds for the minimal and maximal eigenvalues of such selected covariance matrices. Compared with a direct union-bound argument, our results provide substantially sharper guarantees and allow much smaller subset proportions. We further discuss extensions from Gaussian to sub-Gaussian random vectors, and beyond independence to weakly dependent observations, with geometrically strong-mixing Gaussian sequences serving as a representative example of the latter. Finally, we apply the developed concentration inequalities to the K-subspace clustering problem under a low-rank Gaussian mixture model, where the optimal clusters are inherently data-dependent. Our results yield recovery guarantees showing that the clustering error of global minimizers decays polynomially with the signal-to-noise ratio.

We develop a hierarchical Bayesian framework for covariance matrix estimation built on a key observation: while equivariance under the full general linear group GL(p) is well known, it is an extremely restrictive property – estimators equivariant to GL(p) are limited to scalar multiples of the sample covariance matrix and carry considerably larger risks than shrinkage estimators. By contrast, commonly used shrinkage estimators, including the Haff empirical Bayes estimator, and the Ledoit–Wolf estimators, are all equivariant under the smaller orthogonal group O(p). Exploiting this structure, we establish that the Haar measure Bayes rule in an oracle eigenvalue model is the minimum risk estimator within the class of O(p)-equivariant estimators, and derive oracle Bayes rules for the covariance and precision matrices under the squared Frobenius, Stein, and squared Stein loss functions. These oracle rules serve as theoretical benchmarks that dominate all commonly used estimators. To approximate them when the true eigenvalues are unknown, we introduce a hierarchical Bayes model that places a finite P'olya tree prior on the eigenvalue distribution and uses Gibbs sampling to generate posterior draws, yielding both shrinkage estimates for the eigenvalues and approximations to the oracle Bayes rules. Simulations suggest that the finite P'olya tree prior is able to recover the general form of the distribution of the eigenvalues, and confirm that the resulting estimators closely approach oracle performance, substantially outperforming classical competitors for both covariance and precision matrix estimation.

We develop a continuous-time state-space framework for the joint reconstruction of three Cenozoic climate proxies, benthic foraminiferal d18O and d13C and atmospheric CO2, from irregularly and unevenly sampled multi-site, multi-method data spanning the last 67 million years. The latent signals follow a trivariate random walk in continuous time; the measurement equation differentiates the error variance by drill site for the isotopes and by proxy group for CO2, with bias intercepts placing all sources on a common scale, and the transition equation lets the innovation covariance and a deterministic La2004 Milankovitch forcing depend on the prevailing climate state. All parameters are estimated by maximum likelihood through the Kalman filter with diffuse initialization. The estimated cross-proxy correlations reverse sign between the early Cenozoic greenhouse and the icehouse, the orbital sensitivity of the isotopes strengthens as continental ice sheets grow, and the reconstructed CO2 path, reported with calibrated confidence bands, places the atmospheric CO$_2$ thresholds of the major Cenozoic glaciations in relation to present-day concentrations.

We study the problem of model selection among probabilistic forecasting models evaluated on datasets of multiple time series. The performance of a model on a single time series is quantified by the average value (score) of a proper scoring rule over a test set, but extending model selection to data sets of time series requires aggregating these scores. Common approaches either rely on scaling scores and averaging them (mean scaled score) or avoid scaling by using alternative statistics such as mean ranks or win rates. However, these approaches can yield conflicting conclusions. We show that such discrepancies arise from the skewness of the distribution of the scores, which is particularly pronounced when test sets are short. The skewness can cause non-mean criteria (e.g., mean rank, median, win rate) to select misspecified models. In contrast, the mean score is immune from this problem. We further show that, as the size of the test sets increases, all aggregation criteria converge to the same model selection decision, mitigating these discrepancies. Our experiments on intermittent demand time series, including data from the M5 competition, highlight the importance of sufficiently large test sets; the mean scaled score appears to be the more reliable approach, also because empirically we found its decision to remain consistent when different scaling factors are adopted.

Understanding dynamic systems is a central goal in many scientific disciplines. State-space models provide a general framework for studying latent dynamic systems based on indirect observations. However, classical state-space methods require researchers to specify the parametric form of the system dynamics in advance, which can be challenging when the underlying processes are nonlinear and only partially explained by theory. Gaussian process state-space models address this by learning the system dynamics directly from data. However, estimating these models exactly can become computationally infeasible for moderately long time-series. In this paper, we therefore aim to improve the Bayesian estimation of approximate Gaussian process state-space models to make these models more accessible and facilitate the statistical learning of nonlinear dynamic systems in empirical research. To this end, we first propose two modifications to an existing Gibbs sampler for these models that considerably improve its sampling efficiency and convergence. Second, we use a confirmatory factor analysis measurement model, which reduces identifiability issues and allows researchers to impose a specific measurement structure on the model. Third, we provide a systematically validated software implementation of the model and sampler for applied use in empirical research. To validate the sampler, we conducted a simulation-based calibration which showed that the sampler converged reliably across many simulated data sets and produces well-calibrated posterior inferences. We further illustrate how the model can be applied and interpreted using two empirical examples. Together, these contributions provide a practical and validated workflow for learning nonlinear latent dynamics with Gaussian process state-space models.