Papers

We consider multivariate Poisson autoregressive processes evolving on random networks generated by a Stochastic Block Model. This framework extends existing network Poisson autoregressions by incorporating latent community structure and allowing for non-Markovian dependence, both excitatory and inhibitory. We establish a mean-field approximation and prove that the finite system converges to its limit at the optimal rate $N^{-1/2}$ in the quenched setting with overwhelming probability, both on average and in uniform norm. Our analysis yields quantitative error bounds that are explicit in the time horizon and, under suitable stability conditions, remain polynomial or uniform in time, thereby avoiding the exponential growth typically associated with Grönwall-type arguments. As a consequence, the mean-field approximation is shown to hold uniformly over any time horizon. We further demonstrate how the approximation can be exploited to reduce the computational complexity of parametric inference from $\mathcal{O}(N^2T^2)$ to $\mathcal{O}(NT^2)$. The results provide a theoretical foundation for the analysis and inference of count-valued processes on large random networks with community structure.

Multi-step forecasting systems, including modern artificial intelligence (AI) models, often issue overlapping predictions for the same future time from different forecast origins. These predictions can disagree even when produced by the same model, a phenomenon we call cross-origin incoherence. We propose a Bayesian reconciliation framework that treats overlapping forecasts as noisy and potentially biased measurements of a common latent trajectory. The method operates as a model agnostic post-processing layer, requiring no access to the forecasting models architecture, parameters, or training data. A Bayesian state-space model jointly reconciles forecasts targeting the same future time while accounting for horizon-specific bias and uncertainty. Controlled simulations demonstrate substantial reductions in cross-origin incoherence and accurate latent trajectory recovery under the assumed data generating process. An empirical application to daily S and P 500 volatility forecasts reduces 90% prediction interval coverage from 94.8% to 88.3%, closer to nominal coverage, but increases point forecast error. An additional LSTM experiment demonstrates applicability to modern neural forecasting systems.

Generative foundation models are attracting interest for their ability to produce desired outputs from demonstrations given at inference time, without updating parameters. However, since a few demonstrations cannot uniquely identify the intended task, the challenge is how to learn and sample from an output distribution that reflects this task uncertainty. In this work, we theoretically analyze how a Diffusion Transformer (DiT), pretrained across diverse tasks, learns and generates predictive distributions for a new query from demonstrations. We first show that the natural target to generate from finite demonstrations is not an output derived from estimating a single task, but rather a predictive distribution that captures the task uncertainty remaining after observing the demonstrations. We then prove that a DiT can learn this predictive distribution through score estimation, using attention to aggregate information from demonstrations and diffusion to generate samples. Owing to this property, with sufficient pretraining resources and diffusion sampling steps, the resulting DiT achieves the minimax optimal rate over a Hölder class of test-time tasks. These results imply that DiT acts as a statistically grounded in-context generator capable of generating distributions adapted to new tasks while retaining the uncertainty inherent in finite demonstrations.

Approximate Message Passing (AMP) algorithms are attractive as they are computationally efficient and simultaneously admit a precise characterization in terms of the low-dimensional "state-evolution" recursion. In this work, we establish quantitative universality for AMP with centered rank-one sensing matrices $Z_i=(x_ix_i^\top-I_d)/\sqrt d$, where $x_i$ are independent standard Gaussian vectors. These matrices arise in quadratic regression and learning quadratic neural networks. Their normalized vectorizations have isotropic covariance in dimension $p=d(d+1)/2$, but strongly dependent coordinates. For linear observations with Gaussian noise and prescribed smooth, bounded spectral denoisers, we compare rank-one sensing with a covariance-matched isotropic Gaussian ensemble. At a fixed iteration horizon, with $n\leq Cp$ and uniform bounds on the normalized signal energy, initialization, coefficients, and first four denoiser derivatives, normalized signal overlaps, cross-time overlaps, and smooth linear spectral statistics differ by $O(d^{-1/2})$ in every fixed $L^m$. This yields universality of the normalized mean-squared error, fixed spectral moments, and empirical spectral distributions. Whenever the Gaussian observables admit a joint state-evolution limit, the same predictions hold under rank-one sensing. The proof uses columnwise Lindeberg replacement, leave-one-out analysis, smooth spectral expansions, and high-moment concentration. For prescribed smooth spectral denoisers, the result establishes the AMP universality conjectured by Maillard et al. (arXiv:2408.03733).

The mean travel time to work is one of the most cited statistics of the American Community Survey (ACS), yet its distribution, its dependence on mode and departure time and its relation to other measurement bases are rarely set out together. Using ACS 1-year tables, we report that US workers aged 16 and over who did not work from home averaged 27.2 minutes one way in 2024, against 25.0 minutes in 2006, a peak of 27.6 in 2019 and 25.6 in 2021. Dividing published aggregate minutes by workers, commuter rail, long-distance train and ferry riders averaged 71.0 minutes, 2.68 times the 26.5 minutes of solo drivers, and all public transportation 49.5 minutes. In 2024, 9.3% of commuters, about 13.3 million people, travelled 60 minutes or more each way and 3.0% 90 minutes or more; commuters leaving between 5:00 and 5:29 a.m. averaged 35.1 minutes. State means ranged from 33.2 minutes in New York to 17.5 in North Dakota, and big-metro means from 36.5 minutes in New York to 21.5 in Buffalo. Working from home rose from 5.7% of workers in 2019 to 17.9% in 2021 and stood at 13.3% in 2024. The 2022 National Household Travel Survey recorded an average commute trip of 13.43 miles at 25.51 mph. On different definitions, the 2019 EU labour force survey module gave an EU average of 25 minutes and 30 minutes for the United Kingdom, and England's travel diary a commuting trip of 31.5 minutes in 2025 against 26.5 in 2002. The three bases are reported side by side and never pooled. All inputs are public.

We present a unified, transformation-free linear regression framework tailored for ordinal compositional data, such as distributions across educational levels or aggregate Likert-scale survey responses. Traditional log-ratio approaches often obscure interpretability, ignore the inherent ordering of categories, and struggle with boundary zero values. To overcome these limitations, we propose a deterministic model constrained to the space of column-stochastic transformation matrices, naturally preserving the fundamental geometry of the simplex. By adopting the weighted 1-Wasserstein distance as the loss function, our method directly embeds the ordinal nature of the data into the estimation process. We provide a computationally efficient and globally optimal solution explicitly formulated as a Linear Programming (LP) problem. The framework systematically addresses four distinct regression scenarios, introducing a novel ordinal tensor product to rigorously handle interactions across multiple compositional predictors. Furthermore, we equip the model with a suitable regularization strategy and novel diagnostic tools, including a Wasserstein-based coefficient of determination ($R^2_W$) and an Order Preservation Index (OPI). Extensive simulation studies and two empirical applications demonstrate the robust predictive performance and interpretability of the proposed methodology.

Patient world models used for clinical trial simulation can agree on transition kernels and arm-specific risks, yet disagree on the fraction of patients harmed by switching treatment—the counterfactual quantity that matters for intervention-aware reasoning. We audit this reliability gap in a two-stage shared-response SCM: a categorical intermediate health state is followed by common terminal care. Under independent stages, the sharp harm interval has closed-form endpoints for at most three intermediate states, with an exactness boundary at four states. Declared dependence and response-mismatch budgets yield calibrated outer bounds when stage independence or complete mediation is relaxed; in a symmetric three-state model the entire sensitivity frontier is sharp, $[0,\min\{1/2,1/3+(\rho+\delta)/2\}]$, and shows exactly how budgets erase the gain over endpoint-only bounds. Two eight-variable response LPs propagate interventional uncertainty for finite-sample audits. Exact witnesses verify attainability. On public clinical simulators (EpiCare; sepsis), native configurations show little resolved stage dependence and no additional joint-compatibility gain over pairwise transport—honest negative results for reliability claims. All experiments are locally reproducible; guarantees remain conditional on the stated causal model. The results provide a concrete protocol for deciding when a patient world model is safe to trust for counterfactual harm.

Risk regulation imposes directional constraints on scores; we adopt their strict per-input form – the score monotone non-decreasing in every exposure input – as a normative commitment. Deployed pipelines – monotone hand-crafted aggregates feeding sign-constrained gradient boosting – already satisfy it by composition, so constrained-versus-unconstrained comparisons price a guarantee the incumbent has for free. We instead hold admissibility fixed on both sides and measure what learning the aggregation is worth. Our instrument is a recurrent network whose state is classical risk statistics (an exponentially weighted moving average and a high-water mark with learned transforms), monotone by construction in every input and per MC-dropout sample. The central finding, by functional regression, is a subsumption boundary: a learned monotone channel reproduces the geometrically weighted separable family of hand-crafted statistics, one channel per member, to Spearman $\rho \ge 0.996$, approximates window statistics with measurable ceilings, and fails at consecutivity ($\rho = 0.924$) and time localization (0.628), both structural, and at the exposure floor (0.829), a learnability boundary. One explicit admissible basis repairs each failure (rank correlation 1.000). In or near the separable family, learned and engineered aggregation are substitutes, and the learned channel is never statistically behind at full sample size and specified capacity. Its advantages are incumbent-specific: a committed grid pays up to 0.019 AUC in decay regions it leaves uncovered (the learned channel stays within 0.004 of the strongest engineered consumer at every swept point); the highest-dimensional comparator degrades fastest with scarce data; and beyond the training support, grid-fed tree-ensemble scores go flat while a strictly increasing head keeps ranking. No single incumbent is dominated on all three axes.

Two-sample Mendelian randomization with summary statistics has become a popular MR design in Epidemiology. Its ratio estimate has an average causal effect interpretation only under strong untestable assumptions including linearity. We derive what the per-variant ratio, and the pooled and robust estimators built from it, estimate when a continuous exposure has a nonlinear effect, under two conditions on how genetic variation and unmeasured factors shape exposure: no unmeasured common effect modifier, and an exposure control function. Under either condition, each variant's ratio is a weighted average of the slope of the dose-response curve, with weights that depend on how that variant moves the distribution of exposure and that are nonnegative when the exposure distributions at different genotypes do not cross. Variants that shift everyone's exposure equally target nearly the same quantity; variants acting at particular doses or on the spread of exposure target different ones, some outside the range of any individual effect. Pooled estimates therefore have no fixed interpretation, heterogeneity tests and MR-Egger misread nonlinearity as pleiotropy, and summary statistics cannot reveal the problem, even when the dose-response curve is identifiable from individual-level data. For continuous exposures, two-sample estimates rest on a linearity assumption that should be stated and defended; individual-level exposure data resolve the problem and should be prioritized.

The repetition count that works best for a small language model may not remain best at a larger scale. We study this effect in pretraining with a finite target corpus mixed with generic data at a fixed target fraction. On Wikipedia-derived data and Proof-Pile-2, the ranking of measured repetition counts changes with model size, and a 520M Proof-Pile-2 experiment confirms that reducing repetition from sixteen to eight improves loss while using fewer training tokens. We use loss curves from several smaller models to retain a short list of promising repetition counts for evaluation at a larger scale. On PubMed and Caselaw, candidate sets fixed before target-model training retain the lowest-loss measured count on the original evaluation grids at both 200M and 520M. This supports candidate retention as a practical alternative to exact point prediction. We also relate the pruning regression to an empirical scaling model with two opposing repetition-dependent loss terms. A first-order expansion in log model size yields the linear form used by the selection rule, providing a scaling-based interpretation of the candidate-selection procedure.

Pointwise Mutual Information (PMI) is a core measure of testing word association, and Skip-gram with Negative Sampling (SGNS) is essentially a method that implicitly factorizes a shifted PMI matrix. However, a systematic and well-rounded characterization of finite-sample uncertainty in PMI estimation remains absent and imperative to venture into. We provide a statistical framework for PMI estimation and its connection to SGNS. We prove consistency, asymptotic unbiasedness, and asymptotic normality of the empirical PMI estimator, derive its variance via the Delta method, and, applying stochastic approximation theory, obtain a variance decomposition for SGNS-based PMI estimation that separates data variance from optimization variance. Simulation experiments validate the Delta method approximation. Real-data experiments on the Brown Corpus (d = 100) reveal that SGNS systematically deviates from the theoretical relationship PMI + log K. The empirical relationship shows an attenuated PMI coefficient, an amplified log K effect, and a positive intercept, indicating systematic bias. Word analogy validation confirms the models are effective. The failure to validate the variance decomposition under low-dimensional conditions does not diminish its theoretical value; rather, it identifies the unbiasedness assumption as the key bottleneck and clarifies the gap between asymptotic theory and practice, providing implications for both practice and theory.

Objectives: The study designs a simultaneous speed prediction model based on Bayesian Markov Chain Monte Carlo (MCMC) methods for heterogeneous, non-lane-based traffic on urban primary roads, with the aim of determining the effects of vehicle-specific traffic densities on the operating speeds of the six vehicle types. Methods: The traffic data were obtained from videos recorded by unmanned aerial vehicles on five main urban roads in Dhaka, Bangladesh, and were then processed by a deep learning approach in order to derive information on vehicle trajectories, traffic volumes, densities, and space-mean speeds. Six simultaneous equations for speed prediction were developed using a Bayesian MCMC framework, covering standard cars, utility vehicles, heavy vehicles, three-wheelers, two-wheelers, and non-motorized vehicles. The performance of the models was assessed by means of goodness-of-fit measures. Findings: The coefficients of determination for the models varied from 0.922 to 0.956 and the root mean square errors were between 0.50 and 1.56 km/h for the six types of vehicle. Operating speeds decreased steadily as traffic volume increased, although traffic composition had a significant effect on speed when the demand was moderate or high. Operating speeds increased and congestion was delayed when non-motorized vehicles were excluded. The transferability study showed that the models had a reliable ability to make predictions on separate road sections.

In many scientific studies, evaluating the association between two responses requires adjusting for shared covariates. The expected conditional covariance measures this adjusted association and is used in conditional independence testing and treatment effect estimation. Researchers often estimate it using fitted conditional means, but removing the bias caused by the estimation errors of these conditional means remains a difficult task. This difficulty can increase when the true covariate density is unknown and restricted only by known upper and lower bounds. In this article, we study how accurately the expected conditional covariance can be estimated under these density bounds. We propose a multiresolution estimator that reduces bias without estimating the density. Theoretical results show that, for bounded responses and rough regression functions, an unknown density slows the optimal worst-case rate by a power of sample size compared with a known density. We establish that the optimal rate under density bounds alone strictly improves upon the previously conjectured polynomial rate by a sub-polynomial factor, resolving a long-standing conjecture. Our estimator attains this rate up to a logarithmic factor. Through simulations and a real data resampling experiment on diamonds, we show that our method can yield smaller mean squared errors than standard uncorrected estimators.

The theoretical understanding of multi-layer neural networks is largely confined to overparameterized settings, which obscure parameter identifiability and incur high sample complexity. Neural tangent kernel (NTK) provides a general theory for wide networks, but does not offer efficient sample-complexity guarantees. Recent feature-learning results go beyond kernel methods for single-neuron, multi-index, and hierarchical targets. However, the analysis is often restricted to shallow or specific architectures and to the overparameterized regime. We break this paradigm to achieve parameter-level recovery of deep target networks, albeit by using active data queries. Specifically, we study $L$-layer polynomial networks with even degree-$k$ monomial activations and nonnegative higher-layer weights. This structure makes the target network input-convex, while the optimization landscape remains highly nonconvex with respect to the parameters. Leveraging input convexity and active queries, we propose ASPIRE (Active SamPling for Iterative Recovery via Eigendirections), a layerwise sampling-based diagonalization algorithm that recovers all network parameters to $\delta$-accuracy with sample complexity $ \widetilde O_{k,L}\left(d^{L^2+O(L)}\delta^{-2e}\right) $ in polynomial time. To our knowledge, this is the first parameter-recovery guarantee for deep target networks whose exponent grows only polynomially with depth, as well as the first justification for the effectiveness of using high-quality data in neural network training, with a remarkably exponential separation.

Observing a prediction error at one node can help correct predictions elsewhere, but the benefit depends on the residual dependence between nodes. A graph suggests where that dependence might occur, yet does not establish its sign or strength. We develop a Bayesian model of residual covariance to determine how revealed errors should update a fixed feature-based predictor. For a fixed set of revealed nodes, we derive an exact identity for the change in expected squared error under linear residual propagation. The identity characterizes the optimal linear update and expresses the excess error of any other update as an exact quadratic. Under covariance uncertainty, the Bayes-optimal fixed linear update uses the posterior mean covariance. We average over four positive-semidefinite covariance families representing no dependence, positive dependence, negative dependence, and dependence whose sign can alternate with graph distance. Across eight graph benchmarks, four show posterior support within the distance-profile family for negative covariance between neighbors and alternating signs with distance. Among the evaluated methods, ours achieves the smallest worst-case recall gap to the best method on each graph. The method improves squared error over residual propagation on five of eight graphs, but performs substantially worse on one dense heterophilous graph. In these experiments, selecting the highest-posterior family performs similarly to model averaging.

Finite discrete count data are ubiquitous in fields such as biomedicine and social sciences and these data frequently exhibit over-dispersion, rendering the binomial distribution inadequate for modeling. Although the beta-binomial distribution mitigates this by assigning a subjective beta prior to the success probability, its reliance on a pre-specified parametric shape may result in biased representations of prior information, failing to capture the complex and heterogeneous characteristics of real-world data. To address these limitations, this paper proposes a novel Bernstein-binomial model. By leveraging the uniform approximation properties of Bernstein distribution, we construct a highly flexible and general prior framework for the success probability in binomial distribition, which dynamically adapts to complex data structures, avoiding the subjectivity and restrictions associated with conventional unimodal priors. A comprehensive theoretical framework for both frequentist and Bayesian inferences is established. Furthermore, the model is extended to a regression setting to account for covariate effects. Extensive simulation studies and real-world data applications confirm that the proposed Bernstein-binomial model significantly outperforms existing models. It offers superior fitting accuracy, robustness, and out-of-sample predictive performance when modeling over-dispersed counts with complex latent prior structures.

Near-perfect classification of small clinical datasets can hinge on preprocessing, data leakage, missing-data handling, and the evaluation sample. We reproduced a published chronic kidney disease (CKD) classification study on the UCI CKD dataset (400 records, 24 predictors) under explicit leakage controls. A stratified 70:30 development-holdout partition was fixed before any learned preprocessing; imputation, encoding, scaling, supervised feature selection, and PCA were fitted only within training folds. Nine classifiers, three feature representations, and three missing-data strategies were compared by 5x3 repeated stratified cross-validation and ranked on training cross-validation only (ties reported); a frozen set of finalists was evaluated once on the holdout (discrimination, calibration, Brier score), treated as a locked historical internal sample. Ablation, grouped SHAP values, three further seeds, and post hoc sensitivity analyses examined robustness. The top-ranked selected-feature pipeline (KNN imputation, Extra Trees) reached a repeated cross-validation accuracy of 0.9869 and a holdout accuracy of 0.9833 (exact 95% interval 0.9411-0.9980). PCA (27 components, 95.96% of variance) gave only a marginal, seed-dependent advantage over interpretable original-feature representations, and several pipelines were statistically indistinguishable. Nested resampling of the selection procedure gave an outer accuracy of 0.9845, 0.0095 below the top-ranked pipeline on the same folds; missingness indicators alone yielded a ROC area of 0.85-0.88, and removing renal, urinary, and hematologic fields cut accuracy to 0.89. High internal performance persisted after leakage controls but appears driven largely by clinically proximal laboratory and urinalysis variables and CKD-associated missingness; it is not clinical validation, and external, prospective evaluation remains necessary.

Noisy matrix completion is a fundamental problem in statistical learning and has attracted a substantial amount of interest over the last two decades. In a variety of applications, the missingness pattern is highly informative, yet it has received relatively less attention. Most existing methods either overlook this source of information or oversimplify its generating process, with few attempting to model the association between the data matrix and its informative missingness. In this study, we propose a general modeling framework that jointly models the data matrix and its missingness pattern using a pair of low-rank models coupled by a flexible shared linear structure. This joint low-rank approach incorporates the informative missingness to improve matrix completion and can flexibly capture various associations between the two modes of data. We develop an efficient joint estimation procedure for this framework based on projected gradient descent, and establish local convergence guarantees that unveil its computational and statistical properties. We further demonstrate, through extensive simulation studies and real-world data analysis, that our proposed approach outperforms competing methods, achieving substantial improvements in prediction accuracy on missing entries.

We consider inference for network stochastic differential equations in a small-noise, high-dimensional regime, where the diffusion coefficient vanishes with $\ep\to0$ while the network size and parameter dimension may diverge. We develop a minimum-distance estimation framework based on a deterministic measurement model that may omit nuisance terms in the drift, and derive a non-asymptotic error bound decomposing the estimation error into stochastic fluctuation, nuisance discrepancy, and identifiability terms. For a linearly parametrized interaction model under repeated measurements, we obtain explicit non-asymptotic bounds on the nuisance discrepancy and the identifiability gap in terms of graph connectivity and noise level. Under suitable mixing and balance conditions, this yields consistency and convergence rates for the estimator. We further study graph recovery via adaptive Lasso. The penalized estimator is shown to inherit the rate of the preliminary estimator and to recover the true interaction graph consistently. These results provide a theoretical basis for parameter estimation and edge selection in high-dimensional network SDEs with nuisance drift components.

We propose a stochastic optimization framework in which football performance is modeled as the result of an efficiency-driven production process under uncertainty. The model incorporates economic inputs, a latent efficiency variable evolving according to a stochastic differential equation, and a control variable representing investment in sports development. We derive the associated discounted Hamilton–Jacobi–Bellman (HJB) equation, characterize the optimal investment policy, prove a verification theorem under explicit growth conditions, and show that a bounded investment budget is required for the control problem to be well posed. We relate our efficiency measure to the stochastic frontier and data envelopment analysis (SFA/DEA) traditions used elsewhere in sports economics, and clarify what a dynamic control formulation adds relative to those static approaches. The theoretical results are illustrated with a numerical solution of the HJB equation – validated by a grid-refinement convergence study – and Monte Carlo simulations. A cross-sectional analysis of the 24 AFCON 2025 national teams, based on GDP, population, and FIFA ranking data, estimates country-specific conditional log-efficiency residuals using heteroskedasticity-robust inference. We show that these residuals are mechanically and strongly correlated with FIFA points by construction, quantify this dependence, and interpret the results accordingly. Teams such as Senegal and Morocco consistently outperform their GDP-implied potential, whereas others, including South Africa and Tanzania, tend to underperform. While these patterns are broadly consistent with observed AFCON outcomes, the estimated residual also captures unobserved factors, such as institutional quality, diaspora talent, and demographic structure, that are not included in the model.

In this paper, we study the sample complexities of value and policy learning in finite discounted Markov decision processes (MDPs) under recursive entropic risk preferences with risk parameter \(\beta\neq 0\), assuming access to a generative model of the MDP. We provide a refined analysis of model-based risk-sensitive Q-value iteration (MB-RS-QVI), a plug-in model-based method introduced in prior work, and derive \((\varepsilon,\delta)\)-PAC guarantees for both learning the optimal \(Q\)-value function and an \(\varepsilon\)-optimal policy. Our bounds improve the exponential dependence on the effective horizon \(1/(1-\gamma)\) compared with the best existing guarantees for this setting. In particular, they match the existing lower bounds in their exponential dependence on \(|\beta|/(1-\gamma)\), as well as in \(S\), \(A\), \(\varepsilon\), and \(|\beta|\), up to logarithmic factors. Consequently, our analysis removes the exponential gap between the previously known upper and lower bounds, leaving only a polynomial gap in the effective horizon.

Multivariate functional data arise in many modern manufacturing systems, where multiple sensors record densely sampled process trajectories. Monitoring such data is challenging because nominal variation is strongly correlated across samples, sensors, and time, while faults may appear either as isolated deviations or as structured departures concentrated within a limited number of sensor-specific temporal trajectories. We propose two unsupervised sparse tensor decomposition methods that preserve this multimode structure. Entrywise Sparse CP Decomposition (ES-CP) uses an entrywise \(\ell_1\) penalty to identify localized anomalies, whereas Fiberwise Sparse-Group Lasso CP Decomposition (FG-Lasso) combines entrywise and fiberwise penalties to detect both localized deviations and anomalies concentrated within temporal fibers. Both methods represent nominal process behavior through a low-rank CP decomposition and are estimated using alternating optimization with closed-form sparse-component updates. Two simulation studies evaluate performance under different fault structures, signal severities, noise levels, and missing observations. FG-Lasso attains or ties the highest macro F$_1$ score in almost all settings in the first study and achieves the highest macro F$_1$ score. In a multichannel forging-process case study, FG-Lasso and ES-CP obtain macro F$_1$ scores of 0.85 and 0.82, respectively, compared with 0.69 or lower for TRPCA and PCA-based anomaly detectors. The results demonstrate that explicitly matching the sparse penalty to the anticipated fault structure improves both anomaly detection and fault localization in high-dimensional functional processes.

Rugby Union's continuous play, variable-length possessions, and multiple scoring outcomes make it difficult to model with standard approaches. We present a Markov Chain framework that represents Rugby Union as a finite state space over score differential, remaining time, and pitch location, with transition probabilities estimated from historical match data. The resulting model supports simulation and analysis of in-game trajectories, including win probability and expected point estimates from any given state. We validate the model through calibration analysis and comparison to models that predict only final outcomes. Our state-space approach produces well-calibrated, competitive estimates while additionally enabling analysis of how games evolve between states.

We introduce an algorithm for efficiently preconditioning log-concave and log-smooth distributions that scales logarithmically with the condition number of the underlying distribution. Based on Gaussian cooling, our multistage method approximately samples from a sequence of well-conditioned distributions to construct a preconditioner for each subsequent stage of the cooling schedule. This method is motivated by existing practical approaches for mass matrix estimation in Markov chain Monte Carlo (MCMC), in which an effective preconditioner for the target distribution is estimated from samples produced during a warm-up period. However, our approach is conceptually distinct from existing approaches and moreover permits theoretical analysis. Our algorithmic framework is flexible, allowing essentially any log-concave sampling algorithm to act as the subroutine within each cooling stage. For instance, using first-order rejection sampling (FORS) as the sampler, the first-order query complexity of our method is $\tilde{O}(d^{4/3} \log \kappa)$ for a large class of log-concave and log-smooth distributions with condition number $\kappa$ and $\tilde{O}(d^{1/3} \log \kappa)$ for translation-invariant distributions. In the translation-invariant setting, we provide numerical experiments for the lattice $\phi^4$ model, a simple lattice field theory.

We present a Bayesian framework for assessing the calibration of nested uncertainty sets from simulation studies. Calibration is typically evaluated by repeatedly simulating data from a known model and comparing empirical coverage probabilities with their nominal values, an approach that often requires hundreds or thousands of independent model evaluations to obtain reliable estimates. Such computational demands can be prohibitive when the forward models are expensive to evaluate, for example in data assimilation using a large Earth-system model, or may not be applicable, for exmaple when the output of the simulation is a single or multiple uncertainty sets rather than a full distribution. We present a method that leverages the nested structure of uncertainty sets across nominal coverage levels to assess the calibration of these models. The method, which is based on a Dirichlet-Multinomial model, allows for testing both exact and region-based notions of model calibration in an efficient manner using Bayes factors. Unlike methods that require a full predictive distribution, this method works for procedures that directly produce prediction or confidence intervals. We demonstrate the method on a data assimilation problem containing model form error and show that it is capable of detecting miscalibration even when there are only a limited number of simulation studies available.

Sampling from a probability distribution with access only to its log-density is a foundational problem in statistics, machine learning, and the physical sciences. This problem is studied in the data-free score-based diffusion literature, which samples the target density via estimating scores of densities produced by convolution with Gaussian random variables and implementing a reverse SDE. In this framework, score estimators are constructed via Markov chain Monte Carlo, importance sampling, or rejection sampling. The estimators often operate far beyond the regime where mixing time and error control can be established and have large bias and variance when estimating the underlying score. In our method we propose post-processing: splitting a fixed computation budget across multiple evaluations, querying the estimator at different points, and fitting a single pooled estimate via linear combination to minimize bias and variance at a target location. The final value is an affine map of a vector of multiple score estimates. We show this reduces score estimation error when compared to the computation-matched base estimator queried at a single point. We demonstrate improved sample quality over other score estimation methods matched for computation cost on difficult target densities such as log-concave densities with large condition number, multimodal densities, and high-dimensional densities.

Self-evolving large language model (LLM) systems repeatedly propose, evaluate, and incorporate updates to prompts, skills, or other persistent artifacts. Despite their growing effectiveness, these systems typically operate under a predetermined iteration or compute budget, without a principled criterion to determine when further evolution is no longer worthwhile. This can lead to two undesirable consequences: unnecessary computation after performance has saturated and the risk of returning late updates that overfit or exploit the evaluation signal. These issues motivate us to study two fundamental questions: when should a self-evolving system stop, and what should it output once it stops? We address the first by formulating an online sequential testing problem and constructing an anytime-valid restart detector using the per-item paired evaluation outcomes already produced by self-evolving LLM systems. We address the second by formulating a change-point estimation problem and using the estimated transition to select an earlier artifact for output. The resulting procedure is plug-and-play and requires no modification of the underlying self-evolving algorithms. Across two self-evolving frameworks, three LLM model families, and five benchmarks, our method substantially reduces computation costs while maintaining comparable unseen-test performance. For example, on SearchQA with SkillOpt and DeepSeek V4 Flash, our method stops at round 4 rather than the full budget of 40, reducing token usage by 91.6% while achieving 82.43% unseen-test accuracy versus 82.00% under the full-budget run.

We study model-free online reinforcement learning (RL) for episodic tabular Markov decision processes, focusing on both gap-dependent regret and policy switching cost. While fine-grained gap-dependent analysis has been established for model-free RL algorithms using Hoeffding-type exploration bonuses, such results for model-free algorithms with variance-based exploration bonuses remain unknown, despite their superior worst-case and coarse-grained gap-dependent guarantees. In this paper, we resolve this open problem by establishing the first fine-grained gap-dependent regret upper bound for UCB-Bernstein+, a refined UCB-Bernstein algorithm, in variance-aware model-free online RL. Moreover, by integrating a stage-wise policy update design into our fine-grained framework and using refined variance-based bonuses, we achieve the best-known gap-dependent local switching cost to date. In addition, our analysis yields improved worst-case guarantees for both regret and local switching cost over the original UCB-Bernstein algorithm. Numerical experiments further demonstrate that UCB-Bernstein+ achieves favorable empirical performance in both regret and local switching cost.

The log predictive density score (LPDS) and related estimators of the expected log predictive density (ELPD) have become the de facto standard for Bayesian model comparison. The LPDS is a natural out-of-sample measure of predictive density performance that accounts for parameter uncertainty, is less dependent on the prior than the marginal likelihood, and is easily computed sequentially using standard Monte Carlo from posterior draws. Our article demonstrates, however, that the standard Monte Carlo estimate of the LPDS can be extremely variable and severely biased when the test data has outliers, leading to misleading conclusions and suboptimal model choices. Estimating the LPDS as the difference between the log marginal likelihood of the full dataset and the training data is shown to be a dramatically more efficient approach. This estimator is also shown to be a much more effective fallback for observations that have been flagged as problematic by the widely used PSIS-LOO model comparison framework in Stan.

Simulation is a fundamental computational primitive in reinforcement learning (RL), yet conventional simulation explicitly generates individual trajectories even when downstream procedures use only aggregate statistics. To address this, we develop an exact fast-simulation framework for finite-horizon tabular Markov decision processes. Our framework has two complementary modes. In direct batch simulation, a batch is represented by its aggregate Markov flow. With sufficient parallel simulation resources, this flow can be obtained by trajectory aggregation; when such simulation is unavailable or costly but the initial state and transition distributions are directly accessible, we instead generate an identically distributed flow through forward Markov-flow sampling without materializing individual trajectories. The latter reduces the simulator-side computational dependence on batch size $m$ from $O(m)$ to $O(1)$. In adaptive batch simulation, when batch length is determined by a data-dependent condition, exact multivariate-hypergeometric splitting recursively refines a candidate Markov flow while preserving the conditional law, reducing the cost dependence on $m$ from $O(m)$ to $O(\log m)$. Together, these modes accelerate simulation by keeping trajectories aggregated whenever possible and refining flows only when required to locate data-dependent boundaries. The framework applies broadly across simulator-based, offline, and online batch or stage-based RL, as illustrated with representative algorithms from each setting.