Papers

Computational models in science and engineering are often assessed by checking whether the residual norm is consistent with the assumed noise level. This can be misleading in smoothing inverse problems: structured model errors may be attenuated in observation space, leaving residual magnitudes below practitioner discrepancy thresholds while coherent residual patterns remain. As a result, residual-norm diagnostics can accept fitted models that still give biased parameters, predictions, or quantities of interest. We propose a structure-sensitive sequential diagnostic based on e-processes. The method uses a portfolio of spatial residual-pattern experts, updates their likelihood-ratio wealth as observations are processed, and rejects the fitted model when the aggregate wealth crosses a prescribed threshold, giving anytime-valid type-I error control for a fixed fitted model. We compare the method with Morozov discrepancy checks, fixed-sample residual tests, and batch projection tests. Across three inverse problems (elliptic diffusion, two-dimensional Stokes flow, and a glaciological ice-stream inversion implemented in the community finite-element model icepack) we demonstrate how standard discrepancy checks accept misspecified fits that produce materially wrong quantities of interest. Structure-sensitive batch tests detect these failures using the full dataset, while the e-process detects them earlier from a fraction of the observations. After rejection, the expert wealth attributes the evidence to residual patterns in the chosen dictionary and provides a basis for exploratory model correction.

Granular Instrumental Variables (GIV), introduced by Gabaix and Koijen (2024), use idiosyncratic shocks to large units as instruments for aggregate endogenous variables. I develop asymptotic theory for GIV in large panels where both the cross section and time series grow (\(N,T\to\infty\)). I make two main contributions. First, I show that estimating the individual shocks used to construct the granular instrument has a first-order effect on the asymptotic variance of the GIV estimator. Existing methods do not account for this effect, leading to incorrect standard errors. I derive corrected inference. Second, I characterize how instrument strength changes as \(N\) grows, yielding strong, nearly-weak, and weak regimes. In the newly characterized nearly-weak regime, the instrument vanishes asymptotically but the estimator remains consistent and asymptotically normal, at a slower rate, when \(T\) grows sufficiently fast relative to \(N\). I apply the theory to estimate short-run supply and demand elasticities for crude oil, copper, and natural gas.

Grouped orthogonal arrays were introduced to address experimental design problems arising in computer experiments with grouped inputs, as well as in physical experiments where interactions between factors from different groups are assumed to be negligible. Motivated by the growing need for flexible and efficient designs under such settings, this article develops several constructions to expand the existing catalogs of grouped orthogonal arrays. The proposed constructions provide a large collection of new grouped orthogonal arrays with significantly larger numbers of groups and group sizes.

A random geometric graph (RGG) is generated by first sampling latent points $x_1,\ldots,x_n$ independently and uniformly from the unit sphere in $\mathbb{R}^d$, and then connecting each pair $(i,j)$ if $\langle x_i,x_j\rangle$ exceeds some threshold $\tau$. We study the sharp detection threshold – the largest dimension at which the RGG can be statistically distinguished from the Erdős–Rényi graph with the same edge density $p$. This threshold is conjectured to be $d \asymp (nh(p))^3$, where $h(p)=p \log \frac{1}{p} + (1-p) \log \frac{1}{1-p}$ is the binary entropy function. Previous works proved this conjecture for dense graphs with constant $p$ and, up to polylogarithmic factors, very sparse graphs with $p=\Theta(1/n)$. In this paper, we prove that detection is impossible when $d\gg (nh(p))^3$ and $d\ge (1+\epsilon) n$ for any constant $\epsilon>0$, thereby resolving the conjecture in the regime $p\gtrsim n^{-2/3}/\log n$ and improving upon the state of the art in the regime $1/n \ll p \ll n^{-2/3}/\log n$. The key to our proof is a sharp analysis of the posterior distribution of the latent points given the observed graph, obtained through an information-theoretic comparison argument combined with strong log-concavity.

Although neural networks are remarkably effective, their underlying optimization principles remain theoretically elusive, often characterized by non-convex landscapes and stochastic heuristics. In this work, we propose a paradigm shift by replacing the discrete training problem of shallow neural networks with a well-posed continuum variational surrogate. We identify a family of $\lambda$-convex functionals over parameter densities in weighted Sobolev spaces and prove that these variational problems are globally well-posed, stable, and exhibit unexpected almost $C^3$ regularity. Unlike existing Wasserstein-based or Mean-Field approaches, which often face limited regularity and discretization challenges, our formulation provides direct access to elliptic regularity and convex analysis. This allows us to prove that the optimal parameter density can be obtained by solving a single linear system, bypassing iterative optimization entirely. We establish explicit generalization error controls at a rate of $1/\alpha$ relative to the regularization parameter, and prove that finite-width networks of size $N$ achieve the continuum optimum at an $O(1/N)$ rate. This perspective bridges the gap between the Neural Tangent Kernel (NTK) and feature-learning regimes, providing a principled framework for understanding over-parameterization through the lens of variational calculus.

Background. Bayesian kernel machine regression (BKMR) is widely used for exposure-mixture analyses with binary outcomes through a probit extension. Because a bkmr fit can complete without providing adequate effective posterior information, simulation studies should separate execution success from MCMC convergence diagnostics. Methods. We evaluated the public bkmr probit workflow using bkmr::SimData() for data generation, bkmr::kmbayes() for model fitting, and posterior for convergence diagnostics. The balanced generator used family = "binomial", hfun = 2, beta.true = 0.5, ind = 1:2, and M = 4. SimData() generated the covariate as X = 3*cos(z1) + 2*rnorm(n). Four chains were initialized with chain-specific randomized starting values generated reproducibly from the fixed initial-value base seed 20260621. These values affected only the initial state of the sampler and did not alter the BKMR model, default priors, or Metropolis-Hastings proposals. Results. Of 431 prespecified tasks, 430 returned fitted objects and one task had a numerical non-completion. Diagnostic adequacy was limited: rank-normalized R-hat <= 1.01 threshold was achieved in 55/431 tasks, bulk-ESS >= 400 in 85/431, tail-ESS >= 400 in 44/431, and both ESS criteria in 44/431. The primary diagnostic criterion, R-hat at or below the 1.01 threshold with both bulk-ESS and tail-ESS >= 400, was met in 30/431 prespecified tasks, corresponding to 30/430 completed fits. Conclusions. Completion of probit BKMR fits in bkmr should not be equated with convergence of the retained MCMC draws. Applied analyses should report the number of chains, warmup and retained iterations, rank-normalized R-hat, bulk-ESS, and tail-ESS rather than rely on a fixed iteration count or on fit completion alone.

Plausbility is a theoretical framework that allows to conduct exact inference in general parametric families. We introduce R-packages plausibility that implements this framework for a wide class of regression models. Plausibility can also be used to test penalized regression models such as estimated by package glmnet. We illustrate the package using a number of R data sets Through a class-based mechanism, the package can be easily extended. We illustrate and discuss computation aspects of the implementation and their impact on real-data analysis.

We propose MOMENT (MOment-Based Mixed-Effects SelectioN and EsTimation), a stage-wise moment-based framework that exploits second-order cross-moment identities to select and estimate the random-effects covariance matrix and fixed-effects coefficients. By inducing sparsity through its diagonal under a positive semidefinite constraint, the random-effects selection problem reduces to a smooth constrained convex optimization problem that can be solved efficiently by projected gradient descent. We further establish finite-sample theoretical guarantees for the proposed procedure, including random-effects selection consistency and fixed-effects selection consistency under joint sub-Weibull errors. Simulation studies show that MOMENT performs competitively overall and can substantially outperform separate univariate analyses when responses are correlated. An application to the hemodialysis dataset demonstrates that the proposed method yields an interpretable and flexible approach for multivariate longitudinal data.

Suitability scoring for outdoor activities (kitesurfing, paragliding, ski touring) maps environmental conditions to a go/no-go verdict via expert-defined curves. These curves conflate two distinct quantities: the intrinsic difficulty of a condition and the skill of the person facing it. We introduce Inverse Suitability, a continuous-item Item Response Theory (IRT) model that identifies both from behavioural outcomes alone. Each outcome is a triple (rider r, condition metric x at site s, binary outcome y); we model P(y=1) = sigma(a (theta_r - delta(x, s))), where theta_r is latent rider skill, delta(x, s) is a latent difficulty function anchored to a physics-derived expert curve as its prior, and a is a discrimination parameter. The formulation is strictly more general than a single suitability curve, which it recovers exactly when skill is integrated out under the population distribution. Parameters are estimated by marginal maximum likelihood with Gauss-Hermite quadrature; identification holds when the rider-by-condition incidence graph is connected, with a documented single-curve fallback otherwise. We validate via synthetic recovery: on a reference cohort (80 riders times 30 outcomes) the model recovers latent skill at r = 0.96, locates the difficulty minimum within 3 units of ground truth, and improves held-out Brier Skill Score by +0.33 over the expert-curve baseline. The recovered difficulty function defines a measurable, site-level construct, an intrinsic difficulty atlas, that existing meteorological observation networks do not capture. All results reproduce from a single command on synthetic data, requiring no proprietary observations.

We propose a model agnostic methodology to measure lag relevance in machine learning forecasting models applied to univariate time series. Particularly, we are working in the context of time series using the frameworks of Ghost variables and Shapley values, together with additive importance measures, to introduce the auto-relevance and partial auto-relevance functions as the lag importance values. Additionally, we propose a novel method to replace absent features in coalition based methods with a one step forecast from the same model. We evaluate these proposals under different simulations and real data cases. This combined framework perspective is particularly suitable for time series. In addition, to show our discoveries we use a pull of models from the seasonal ARMA family and recurrent neural networks. We found that the calculated relevance measures successfully demonstrate the expected lag structure in almost all cases.

The classical $k$-means clustering cannot be directly used to incomplete data, and existing $k$-means-based clustering for missing data primarily focus on improving the practical accuracy of clustering, whereas most of them lack theoretical guarantees in the asymptotic sense. In this paper, we investigate the statistical properties of $k$-means clustering in the presence of missing data. We first establish the $\sqrt{n}$-excess risk bound and prove the consistency of the estimated cluster centers under general missing mechanisms. For the Missing Completely at Random (MCAR) mechanism, we further derive the $\sqrt{n}$-convergence rate and asymptotic normality of the estimated cluster centers. Moreover, we study in what cases the cluster centers estimated by incomplete data converge to the true cluster centers of original fully observed data, and give a sufficient condition about the missing probability and the separation among true clusters. These results provide a theoretical guarantee for missing-data-$k$-means. Notably, our analysis reveal that under MCAR mechanism, both achieving the $\sqrt{n}$-rate and converging to the true cluster centers require $k$ true centers to be distinct in every dimension, highlighting the significant challenges of application in high-dimensional regimes. Finally, we conduct numerical simulations on synthetic incomplete datasets to support our theoretical analysis results.

We show that $L_2$-perturbation theory can be used to transfer rates of convergence in the supremum norm as well as weak convergence in the space of continuous functions from covariance kernel estimators to the associated functional principle components (FPCs). As an application we obtain optimal rates of convergence in sup-norm, including minimax-lower bounds, as well as asymptotic normality for estimating the FPCs in a discrete observational model with errors under fixed, synchronous design. The sparse to dense transition which has previously been observed for mean function and covariance kernel estimators also applies to the FPCs. Surprisingly, eigenvalue estimation exhibits a discretization-dominated regime under sparse designs, too. Our results further apply to estimators of cross-covariance and long-run covariance kernels, as well as to covariance kernels of derivative processes. We also present results of numerical experiments in which we use the Nyström method to compute FPCs and eigenvalues, and give an empirical illustration to series of daily temperature curves.

In Bayesian statistics, the Laplace approximation provides a computationally efficient approximation to posterior distributions. However, its Gaussian form restricts it to elliptical shapes, limiting its ability to capture important posterior features such as skewness, heavy tails, and narrow high-probability regions. Recent work has addressed this limitation by exploiting Riemannian geometry to push forward Gaussian distributions from the tangent space to a chosen manifold, referred to wrapped Gaussians. While offering greater flexibility, they introduce substantial computational challenges. Sampling requires solving geodesic equations through the exponential map and density evaluation additionally depends on the logarithmic map and Jacobi fields, involving costly differential equation solvers alongside full matrices inverse operations, Christoffel symbols and curvature tensors. To overcome these limitations, we employ the theory of contrast functions to derive tractable approximations of the logarithmic and exponential maps on statistical manifolds endowed with the Fisher–Rao metric and the prior distribution geometry. The resulting methodology bypass the need to compute these geometric quantities and numerical solvers thereby removing the principal computational bottlenecks of existing wrapped Gaussian approaches. Empirical results across a range of models demonstrate that the proposed approximation captures complex posterior geometries while remaining orders of magnitude faster than current geometric state-of-the-art approximation.

We study ridge-regularized log-density-ratio estimation in the Gaussian location model with a common covariance matrix. By affine invariance, the model is written as q $\sim$ N(0, I), p $\sim$ N($\Delta$, I), with linear features, where $\Delta$ is a mean vector. The variational estimator is the empirical Kullback-Leibler (KL) log-normalized fit with a squared L2-penalty on its nonconstant coefficient, and the spectral estimator recently introduced in [1] replaces a single variational problem by a continuum of ridge-regularized least-squares problems. We derive high-dimensional deterministic asymptotic equivalents when the numbers of observations and dimension tend to infinity with fixed ratios. The regularized variational limit is characterized by a scalar entropy minimization problem derived from the convex-Gaussian-min-max theorem (CGMT), while the regularized spectral limit follows from deterministic equivalents for resolvents of weighted sums of two independent Gaussian sample covariance matrices. We use these formulas to compare population risks, with experiments focused on fixed-signal aspect-ratio sweeps and optimized regularization. Our conclusion is that with many observations, under the criteria and asymptotic regimes analyzed here, the well-specified variational estimator has the smaller risk, while with fewer observations, the spectral estimator is favored because its covariance-based construction has lower variance. We also study how a nuclear penalty can be used and partially analyzed to perform feature learning.

Blockwise missingness in multimodal data is usually treated as an incomplete-input problem. We instead focus on prediction for a prespecified observed-modality pattern, where the observed modality set determines the information on which the prediction rule can condition. A procedure that imputes missing modalities, zero-fills unobserved modalities, or trains a single pooled predictor may borrow information across patterns, but it can also mix pattern-specific prediction rules. We propose Multimodal Overlap-aware Shared-specific Alignment and Inter-pattern Calibration (MOSAIC), a pattern-calibrated framework for borrowing across missingness patterns without collapsing their prediction rules. MOSAIC learns shared and modality-specific representations, uses the available representations that overlap with the target pattern to fit a first-stage predictor, and then estimates the calibration gap from target-pattern data. We establish non-asymptotic bounds that decompose the error into overlap effective sample size, calibration gap, and representation-learning error, clarifying when cross-pattern borrowing improves over local fitting and when the improvement is controlled by rule mismatch or representation-learning error. Simulations examine representation recovery and target-pattern correction, and applications to ICU mortality prediction, emotion recognition, and glaucoma classification show gains when target-pattern samples are limited or pattern-specific rules differ.

Many representation learning problems involve directed relations, such as lexical entailment, sentence entailment, ontology hierarchy, and citation links. Standard Euclidean, cosine, and Mahalanobis heads are symmetric, while generic neural scorers can model directionality but provide limited geometric structure. This paper proposes a role-aware neural convex divergence head for asymmetric representation learning. The head applies source- and target-role projections before evaluating an input-convex neural Bregman divergence, yielding a nonnegative structured score in the role-projected space. We characterize its projected-space identity, source-role convexity, directional-gap decomposition, and Hessian-based local curvature. Experiments on lexical, sentence, ontology, and directed graph benchmarks compare symmetric distances, unstructured asymmetric scorers, order/hyperbolic baselines, plain ICNN-Bregman heads, and the proposed role-aware variant. Across ten random seeds on the main semantic and ontology benchmarks, role-aware projections consistently improve directional accuracy over plain ICNN-Bregman heads while preserving zero observed negative divergence rate. The results also identify a boundary case: on large fixed-feature citation prediction, specialized symmetric or hyperbolic baselines remain stronger in ranking accuracy. Overall, the proposed head is best understood as a structured and interpretable plug-in distance module for tasks where directional relations matter.

Despite increasing scale and resolution, many biological measurements remain destructive, revealing only spatial information rather than the dynamics it encodes. By combining flexible representations with mechanistic constraints, physics-informed machine learning offers a promising route to inferring these dynamics from static snapshots. Motivated by subcellular imaging of gene expression, we ask when a static spatial pattern of molecules can identify spatially varying diffusivity, creation, destruction, and boundary exchange, and how different inference schemes perform on the task. A structural identifiability analysis shows that distributed sources are non-identifiable, whereas a point source such as a transcription site can restore identifiability. These limits are further shaped by seemingly innocuous modeling choices: the boundary conditions, the spatial regularity of the underlying dynamics, and even the stochastic calculus convention. We then adapt several physics-informed schemes, differing in how they represent the solution and enforce the governing equations, and demonstrate effective inference from a single snapshot. Physics-informed approaches can thus recover spatial heterogeneities of biological dynamics from static data, but their use should be accompanied and guided by careful identifiability analysis for meaningful interpretation of the results.

Reinforcement Learning (RL) is a sequential decision-making framework in which an agent learns optimal policies through interaction with an environment by maximizing cumulative rewards. Among RL methods, Bayesian Reinforcement Learning (BRL) addresses common practical challenges related to data scarcity by leveraging prior knowledge about the environment and sequential belief updates. However, most BRL approaches require an explicit likelihood function, which is frequently inaccessible or intractable in real-world settings. We propose Likelihood-Free Iterated Batch Importance Sampling (LF-IBIS), a novel algorithm for BRL that updates the agent's beliefs online as new interactions become available. By combining Approximate Bayesian Computation with Iterated Batch Importance Sampling, LF-IBIS enables full Bayesian inference in settings where the environment dynamics are not described by an explicit or tractable likelihood. The method yields approximate posterior distributions over both environment parameters and optimal policies, providing a quantification of policy uncertainty useful for a Bayesian treatment of the exploration-exploitation trade-off. We test the method on a simulation study in response-adaptive randomization in clinical trials, where closed-form posteriors enable validation. Additional experiments address settings where the posterior has no closed form and illustrate online policy updating based on the posterior distribution of the optimal policy.

Football match outcome prediction is a challenging problem because team strength changes over time, match outcomes contain a high level of randomness, and draws play a central role in the result structure. Classical rating systems such as Elo provide simple and interpretable dynamic summaries of team ability, but they do not explicitly model uncertainty and often ignore football-specific contextual information. This paper proposes an adaptive Glicko-2-based rating framework for probabilistic football forecasting and leaguelevel season simulation. The proposed framework extends the standard Glicko-2 model by incorporating football-specific mechanisms, including margin-of-victory adjustment, dominance weighting, structural shocks, home advantage modelling, and an ordered-logit draw model. The framework estimates latent team strength dynamically, converts rating differences into win-draw-loss probabilities, and uses these probabilities to simulate the remaining part of a league season through Monte Carlo sampling.

Identifying heterogeneous populations across which exposure effects vary is essential for transportability applications, cost-benefit analyses, and intervention prioritization. Traditional methods for heterogeneity analyses rely on parametric regression with prespecified subgroups, which may fail to capture complex patterns of effect modification. While recent data-adaptive methods improve high-dimensional heterogeneous effect prediction, they add methodological complexity to analyses and may offer limited insight into key drivers of heterogeneity. In this paper, we propose a novel, conceptual approach for heterogeneity analyses that considers how exposure effects would differ in populations with different compositions by modeling the population-level effect surface as a function of the distribution of effect modifiers. The approach consists of three steps: i) selecting confounders and effect modifiers based on prior knowledge (or alternatively using data-adaptive methods to learn effect modifiers), ii) estimating exposure effects in hypothetical populations with different effect modifier prevalences using transportability methods, and iii) modeling the estimated effects as a function of prevalence values. This approach provides two types of outputs: estimation of the change in the population-level exposure effects attributable to increases in effect modifier prevalence and ranking of effect estimates across multiple effect modifiers and prevalences to identify population characteristics most strongly associated with differential vulnerability. We demonstrate the approach using Demographic and Health Surveys data to examine heterogeneous effects of drought on child stunting and provide a Shiny application to implement this approach in any setting.

We introduce a new copula class, called Lancaster copulas, built from orthogonal expansions of continuous Lancaster probabilities. We derive infinite-series representations for the copula and its density, study truncation effects, and show in numerical experiments that low-order truncations already provide accurate approximation.

Multiple long-term conditions (MLTCs) arise through complex, time-dependent interactions among diseases, yet existing methods often struggle to jointly model disease progression, multimorbidity networks, and high-dimensional risk factors. We propose a structured Bayesian continuous-time Bayesian network (CTBN) framework for learning directed disease-dependency networks from longitudinal electronic health records. The model allows disease transition intensities to depend on existing conditions, pairwise disease interactions, and exogenous covariates. To control the combinatorial growth of interaction parameters, we introduce order-dependent shrinkage priors that increasingly penalise higher-order effects while preserving clinically interpretable main effects. We compare four sparsity-inducing priors, spike-and-slab, structured normal, Bayesian LASSO, and regularised horseshoe through extensive simulation studies. Across multiple data-generating scenarios, the spike-and-slab prior achieved the best network recovery, variable-selection accuracy, and false-discovery control, while continuous shrinkage priors were less effective for hard variable selection. The proposed framework was applied to UK Biobank primary care records, focusing on data from 33,558 participants who were free of the ten selected most prevalent conditions at age 40 and who subsequently developed at least one of these conditions during the follow-up period. The selected spike-and-slab model identified two dominant disease modules: a cardiometabolic cluster centred on diabetes and an inflammatory cluster linking respiratory and atopic conditions.

Bootstrap inference is a cornerstone of computational statistics, but classical bootstrap and BCa intervals can be unstable for nonsmooth estimators such as empirical quantiles, sample maxima, eigenvalue ratios and maximum correlations. The m-out-of-n bootstrap can reduce some failure modes, but it is slow and highly sensitive to the chosen subsampling size. This work develops an adaptive Bag of Little Bootstraps (BLB) procedure in which the exponent {\gamma} in m=n^{\gamma}is selected by a data-driven criterion combining a plug-in asymptotic mean squared error (AMSE) principle with a stability-risk proxy for the BLB variance estimator. The Python package robustboot implements the proposed method together with BCa comparators, nonsmooth statistics, tests, and reproducible examples for quantile inference, eigenvalue-ratio inference, and maximum-correlation screening. We provide an explicit algorithm, theoretical assumptions, a consistency result on a finite candidate grid, extensive Monte Carlo validation across multiple sample sizes, and sensitivity diagnostics for the selected exponent. Simulation studies show that the adaptive BLB improves interval coverage relative to the ordinary bootstrap while preserving computational efficiency. A curated aggregate cholera surveillance illustration from the Democratic Republic of the Congo demonstrates robust uncertainty quantification for public-health thresholds without overstating reconstructed data.

Pretrained generative foundation models cast forecasting as conditional generation from a learned predictive distribution and forecast unseen series zero-shot. We establish three results that turn their reported success into an actionable, mechanistic understanding. First (a positive benchmark result): on a power-controlled study of 1728 rolling-origin forecasts over 36 series from 19 datasets spanning the full range of STL trend strength (F_T in [0.17, 1.00]), a zero-shot Chronos model significantly outperforms four strong classical baselines – drift, seasonal-naive, Theta, and additive Holt-Winters/ETS – with the best mean MASE (1.187 vs. Theta 1.337, ETS 1.656; Friedman chi^2 = 46.08, p = 8.75e-09; Holm-corrected Wilcoxon p <= 0.015 against every baseline; a Nemenyi critical difference separating it from the classical pack). Second (a novel, quantified mechanism): a controlled synthetic experiment with a known trend-generating process shows why – and reveals that the win does not come from better trend extrapolation. When the true trend is linear, damped, or exponential, additive ETS tracks the slope (slope-tracking ratio 1.02, 1.34, 0.98) whereas Chronos systematically under-extrapolates, behaving as a trend-shrinkage estimator (ratio 0.80, 0.49, 0.36). Third (an actionable selection rule): because the advantage is a shrinkage effect, it is predictable from trend strength alone – the generative model wins 78% of low-trend series but only 44% of high-trend ones, and its edge over ETS is significant on the low-trend stratum (0.982 vs. 1.671, p = 0.002) yet a tie on the high-trend stratum (p = 0.18). Trend strength, computable before forecasting from the training context alone, is therefore a practical a-priori indicator of when to deploy a foundation model. We additionally document a calibration shortfall (80% intervals cover 0.77).

Prior work has shown that transformers can perform exact Bayesian filtering within a fixed hypothesis class. Can they also perform Bayesian model selection – identifying the correct hypothesis class from data? We introduce model-selection Bayesian wind tunnels: controlled environments where ground-truth posteriors over hypothesis classes are available in closed form. Using fixed-point-free involutions – whose defining property f(f(x))=x is purely relational – a 2.8M-parameter transformer achieves 0.01-bit entropy agreement with the Bayesian optimum (3 seeds), with both integer tokens and opaque symbols whose meanings change every episode. This extends to non-nested comparisons: involutions vs. 3-cycles (where neither class is a subset of the other) achieve class-posterior MAE under 0.001, demonstrating genuine model selection beyond simplicity/subset bias. We then identify a sharp perceptual access condition: when the discriminative statistic requires arithmetic – modular addition (rotations) or multiplication (f(x)=cx mod p) – model selection succeeds with integer tokens but fails completely with opaque symbols, and this boundary persists under 112x scaling (2.8M to 316M parameters). A stationarity control confirms the operative factor: opaque tokens with a fixed relabeling succeed (0.009-bit MAE), showing that stable semantics, not integer identity, enable circuit compilation. Header subtask diagnostics localize the failure to the composition of header inversion with arithmetic rather than header parsing itself. Probing frontier LLMs on the same tasks shows qualitative Bayesian behavior but a large calibration gap (~55x), measured through lossy probes and therefore directional rather than exact.

Subquadratic alternatives to attention require compromises when applied to multi-dimensional data: standard convolutions lack global receptive fields and input dependency, while recurrent models require rasterizing data such as images, volumes, and partial differential equation (PDE) into an ad-hoc $1\rm D$ scan order that violates their spatial structure. We introduce HyenaND, a subquadratic, global, input-dependent operator that acts directly on the native geometry of multidimensional data through convolutions with implicitly parametrized global, input-dependent multi-dimensional convolutional kernels. Our CUDA implementation, nSubQ, fuses the FFT-convolution path to turn HyenaND's $\mathcal{O}(L \log L)$ scaling into wall-clock speedups. Across long-context genomics, computer vision, medical imaging, and PDE modeling, pure HyenaND stacks match the accuracy of strong attention baselines, while hybrid configurations that interleave HyenaND and attention layers outperform both pure attention and strong recurrence-based hybrids.

The inverse scattering transform (IST) provides the standard theoretical framework for deriving soliton dynamics. Traditionally, such derivations have been of an analytical, rather than data-driven, nature. In this paper, we combine the conceptual framework of the IST with weak-form system identification methods to discover effective soliton dynamics directly from observed scattering data, without assuming prior knowledge of the scattering equations. Our method avoids parameterizing solitary waves via ad hoc curve-fitting by working in the scattering domain, yielding interpretable low-dimensional models that remain valid in perturbed and near-integrable regimes. We demonstrate the performance of the proposed approach on synthetic and experimental data governed by shallow-water equations of Korteweg–de Vries-type and recover models that are consistent with canonical IST theory.

Grid-connected photovoltaic systems (GCPVS) are generally robust but remain susceptible to faults that can compromise energy conversion efficiency or raise safety concerns. Promptly and automatically detecting such anomalies is therefore essential for maintaining system reliability and performance. However, in practice, labeled fault data are rarely available in real-world deployments, which limits the applicability of supervised approaches. Conventional unsupervised baseline models, including a one-class support vector machine (OCSVM), isolation forest (iForest), and local outlier factor (LOF), are trained on normal operation data and assign anomaly scores reflecting how closely new observations resemble that baseline. Although these methods already accommodate non-linear behavior to varying degrees, kernel-based formulations offer further flexibility in shaping the decision boundary; however, tuning the kernel hyperparameters ordinarily requires some prior knowledge of the fault regime. We overcome this limitation by proposing kernel-based norm monitoring (KNM), a non-linear, unsupervised, window-based fault-detection method designed for continuous processes. Although the paper focuses on the GCPVS as a case study, KNM is a general-purpose monitoring framework applicable to a wide range of industrial processes. Using the Grid-connected PV System Faults (GPVS-Faults) dataset operating in intermediate power point tracking (IPPT) mode, KNM is evaluated in two fault scenarios, sensor faults and partial shading, against three benchmark techniques: OCSVM, iForest, and LOF. KNM achieves up to 99.1% and 98.3% accuracy on the two fault scenarios, respectively, using the Cauchy kernel, compared to 93.5% for the best-performing benchmark. The method is interpretable, and variable contribution plots are proposed to support fault identification.

Empirical research rarely admits a unique analysis. Different analytical choices can lead to different conclusions from the same data, yet these hidden forking paths are difficult to observe. We show that AI agents capture much of the analytical variation among human researchers while making these paths explicit. Across four high-stakes domains, assigning different personas is sufficient for AI agents to report divergent, often opposing, conclusions from the same data and question, with findings systematically aligned with those beliefs. In a study in which 42 human research teams analyzed the same immigration dataset, AI agents reproduced 72% of the human ideological gap in reported effect estimates. Despite reaching opposing conclusions, it is difficult to identify clear issues in each analysis based on the final AI reports: 86% passed independent AI review and 78% passed majority human expert review. These findings suggest that the central challenge is often not flawed analyses, but selective exploration and reporting from a large space of methodologically defensible analyses. AI agents may amplify this longstanding problem by making such exploration inexpensive and scalable. To address this, we introduce the m-value (multiverse value), the probability that an analysis path would produce a claim at least as extreme as the reported one. We further introduce Agentic Bootstrap, which estimates the m-value by using AI agents to sample plausible analysis paths. Applied to the human immigration study, 13.5% of reported human analyses fell in the most extreme 5% of the analysis space (m<0.05). Scientific evidence should therefore be evaluated not only by a single reported analysis but also by its position within the distribution of analyses that could reasonably have been reported. Agentic Bootstrap makes this distribution observable and turns it into a criterion for scientific credibility.

Recent work has established a fundamental trilemma between Byzantine robustness, local differential privacy (LDP), and optimization error in distributed learning. We show that this trilemma does not universally extend to generalization error, but instead depends critically on the privacy regime. Specifically, in the high-noise regime (strong privacy), we prove that increasing privacy reduces the generalization error, i.e., there is no tension between robustness and privacy. In the low-noise regime (weaker privacy), however, the tension between robustness and privacy reappears and increasing privacy indeed degrades generalization. Our theory explains this surprising non-monotonic behavior of the generalization error via matching lower and upper bounds on the algorithmic stability of Byzantine-robust distributed learning under LDP constraints. We corroborate and further analyze these theoretical findings with empirical evaluations.