Papers

This paper contains new 80-digit positive-weight Gauss-Hermite and positive-weight interior-node Gauss-Legendre quadrature rules for up to five dimensions and varying polynomial degree accuracy (depending on quadrature type and dimension). Some of these rules improve on the best available rules in the literature and some offer rules where none (other than the tensor product) existed. The results were produced by combining methodology developed by previous researchers with new approaches. The full precision rules themselves are at https://doi.org/10.5281/zenodo.22881864. Software in Julia, Python, and R providing the rules is available via a package registry and/or GitHub: Quadriceps.jl (both 64-bit and 128-bit; https://github.com/NittanyLion/Quadriceps.jl), quadriceps-py (64-bit only; https://github.com/NittanyLion/quadriceps-py), and quadriceps-r (64-bit only; https://github.com/NittanyLion/quadriceps-r). The software used to create these rules is available via PositiveWeightQuadratureSolvers.jl (https://github.com/NittanyLion/PositiveWeightQuadratureSolvers.jl). The replication package is at QuadricepsReplicationPackage.jl (https://github.com/NittanyLion/QuadricepsReplicationPackage.jl). A snapshot of all five packages is archived at https://doi.org/10.5281/zenodo.22883240.

We study decentralized partially observable team decision problems with low-rank latent dynamics and unknown system models. The proposed framework combines team-theoretic equivalence with low-rank model representations to address cooperative decision-making in partially observable Markov decision processes without prior knowledge of the transition model. Each team member makes decisions based on local private information and delayed common information shared across the team. Using only this available information, each member learns an approximate low-rank Markov decision process and applies least-squares value iteration to compute its policy. This yields a fully decentralized learning and planning algorithm that requires neither a centralized coordinator nor centralized training. We show that the resulting member-side solutions approximate the centralized team solution: despite partial observability, unknown dynamics, and delayed common information, each member recovers the corresponding component of an approximate team-optimal policy. We further establish finite-sample performance guarantees and derive a corresponding sample-complexity bound for the proposed algorithm.

Birnbaum's theorem states that the sufficiency and conditionality principles jointly imply, and are implied by, the likelihood principle. To enable formal certification of the results in Lean, we make the conditionality relation explicit through a finite-sample formulation based on \cite{Evans2013}. We reformulate the theorem as a statement about classes of statistical procedures that preserve statistical relations, and show that the proof reduces to the propagation of evidential equality along chains of conditionality-related inference bases (experiment–observation pairs). The main results are illustrated by Venn diagrams, and a worked finite example exhibits two pairs of inference bases: one related by conditionality but not by sufficiency, and the other related by sufficiency but not by conditionality. For finite parameter spaces with at least two points, the class of likelihood-invariant procedures with codomain $[0,1]$ is strictly contained in the class of sufficiency-invariant procedures.

We consider simultaneously testing hypotheses about thousands of units, e.g., genes or proteins, where each unit yields a handful of replicate measurements and we test whether its mean is zero. A widely used approach in genomics, implemented in the limma software, borrows strength across units to learn the distribution of the unit-specific variances, then computes moderated t-statistics. Here we develop the first procedures using moderated t-statistics that (i) control the false discovery rate (FDR) in finite samples under independence across units and null group invariance, without assuming limma's hierarchical model, and (ii) match the power of an oracle local false discovery rate procedure in a sparse asymptotic regime under limma's working model. Benjamini-Hochberg (BH) with standard t-test p-values has asymptotically zero power in the same regime. Our approach learns the variance distribution from statistics that depend on each unit's data only through its orbit under a compact group of transformations that preserves the null distributions, such as sign flips or orthogonal rotations. We then calibrate the resulting moderated t-statistics against group-transformed statistics pooled across units to obtain compound p-values, which we use with BH and a close variant. For small finite groups, we also construct a Selective SeqStep+ procedure using the same learned scores. Our approach extends to two-sample tests and tests of linear model coefficients.

Clustering is an unsupervised learning technique that partitions unlabeled data into groups. Most existing methods require user-specified parameters, such as the number of clusters or neighborhood size. Conversely, we propose automatic depth-based local center clustering (A-DLCC), a fully data-driven method that eliminates numerical parameter tuning. A-DLCC uses the $\beta$-integrated local depth to identify stable exemplars, points consistently central across multiple locality levels, termed local centers, which are ranked by their representativeness. Each local center induces a group of similar points, with group-level similarity measured by a proposed nonparametric metric called group-level local similarity. To guide merging, we incorporate the bottleneck path idea from graph theory, which forms the basis of our adaptive merging criterion. Based on this criterion, we design a single agglomeration rule in which a group is either absorbed by a neighbor it reaches better than itself or bonded to a neighbor that both sides find more reachable than their own background, every merge being additionally required to be carried by a contact stronger than a configuration-model null expects. The rule automatically estimates the number of clusters and decides when to stop merging. Experiments on synthetic and real data show that A-DLCC produces interpretable clustering results without parameter tuning.

We develop an exact simulation algorithm for scalar diffusion paths and diffusion bridges when the Poisson potential is unbounded in both tails. The method reconstructs the realized range of a Brownian bridge proposal by sampling its maximum and location, together with the maxima and locations of the two adjacent restricted Brownian meanders. Conditional on this finite information, the remaining path decomposes into four conditionally independent interval-constrained Brownian bridges, which can be sampled exactly at the Poisson times required by the rejection test. In contrast to constructions based on an enclosing range layer, the proposed representation retains the exact extrema and their locations. Our algorithm returns an exact finite-dimensional skeleton without time-discretization error and permits exact post-acceptance refinement at arbitrary finite collections of times. Numerical experiments validate the resulting finite-dimensional laws and identify the restricted-meander extremum simulation as the principal computational cost in the nonlinear example.

With the rapid advancement in data collection devices and storage capacities, we have access to increasing amount of spatial point pattern data where each event location is augmented by multiple, potentially non-scalar marks. Therefore, there is need for efficient analysis techniques to investigate the structural relationships between the marks. In this paper, we recall the distance covariance and distance correlation and adjust them to the marked point process setting. As a result, we introduce a novel class of mark characteristics for single real-valued marks as well as multivariate combinations of marks, including mixtures of integer- and real-valued quantities, and non-scalar marks.

Offline reinforcement learning and off-policy evaluation evaluates dynamic treatment rules based on retrospectively collected data prior to deployment. In recent AI applications, state and reward information is recorded as complex text or image, which recent AI advancements such as LLM-as-a-judge can label with unknown bias. Expert annotation may be available but at a higher cost. For example, safety classification via cheap but imperfect classifiers vs. expensive expert review. We show how a limited budget for ground-truth data-annotation can be used via doubly-robust OPE with missing rewards, and we optimize variance-optimal annotation probabilities for sequential off-policy evaluation, where the target policy value is estimated from annotated data. We characterize the optimal annotation probabilities for sequential forward-monotone annotation protocols, and provide a feasible batch-adaptive implementation. Our work is motivated by a collaboration with a homelessness services nonprofit that writes casenotes for individuals over time. Our method can be used to unlock trustworthy inference from casenote data and answer new inferential questions such as: how does expanding outreach effort over time affect progress towards a housing application and improvement in housing placement? In simulations and on two real datasets - casenotes from the nonprofit and human-preference votes from LMArena - we see reductions in RMSE of 34-65% for housing placement and 17-68% for progress towards a housing application at budgets of 40% of full annotation and above, and by 55-62% at every budget on LMArena.

Commonly, classifiers and monitoring procedures are trained from labeled data by optimizing an objective such as the misclassification rate. This may lead to unrepresentative conditional distributions of the outcome (the labels) given important external variables, different from the conditional laws in the population. We show how to modify any given threshold-type classifier resp. monitoring rule to achieve representative conditional label prediction by using adapting the threshold to a covariate $Z$ (the context) to distribute sensitivity while maintaining the false alarm rate. In case that the alarm event is unknown, this approach also allows to (approximately) infer the event in terms of a thresholding rule. The approach is implemented by a computationally cheap nonparametric estimation procedure, and its properties are studied in terms of nonasymptotic error bounds and asymptotic distribution theory including empirical process theory. These results allow to construct uniform confidence bands, functional hypothesis tests and change-detection procedures. For the well known FICOS credit scoring example, often used in interpretable machine learning, threshold adaptation leads to an easily interpretable decision rule which can compete with state of the art methods including transformers, in terms of common classification metrics.

This paper studies adaptive experimental design for multiple testing, where an experimenter sequentially chooses which hypothesis to sample. We propose the e-value-based posterior sampling (e-PS) procedure, which uses the empirical average of log e-value increments to guide randomized sampling and applies e-BH to construct rejection sets. Under conditionally valid e-value increments, the procedure controls the false discovery rate at arbitrary stopping times and produces nested rejection sets. We establish high-probability bounds on the number of samples needed to discover all nonnull hypotheses in terms of the growth and concentration of the underlying e-processes. We specialize these bounds to simple-versus-simple, composite-versus-simple, and simple-versus-composite testing. Simulations and experiments using joke ratings and watermarked text illustrate the procedure's power under limited sampling budgets.

We study statistical rates in entropic optimal transport in the semi-discrete regime where one measure has finite support and the other is subGaussian. Our main result establishes parametric convergence rates for the empirical dual potentials to their population counterparts, with no dimension dependence in the leading term. Our result relies on tailored strong concavity analysis of the semi-dual objective, coupled with specialized bounds for the semi-discrete potentials. As a consequence, we obtain fast rates for downstream quantities derived from the optimal coupling. Chiefly, the empirical barycentric projection achieves a squared-error rate $n^{-1}$, matching the fully compact case and improving over the less favorable $n^{-1/2}$ rate known for fully subGaussian settings. Altogether, these results may indicate a lower complexity adaptation phenomenon whereby the statistical complexity of the barycentric projection is governed by the discrete measure. As an application, we analyze Sinkhorn-EM, an EM-type algorithm in which the E-step is replaced by an entropic optimal transport problem. In a well-specified and balanced two-component Gaussian mixture model, we prove $\sqrt{n}$-consistency of the empirical iterates to their population counterparts for any fixed number of iterations, matching classical EM rates up to a $\sqrt{\log n}$ factor. Simulations support the theory.

Sparse Gaussian processes achieve $O(N)$ inference by replacing the kernel with an appropriate expansion in a fixed basis $\{\phi_j\}$ on the input space. Given a compute budget $M \ll N$, practitioners conventionally truncate the basis to its first $M$ entries. Nothing in the formalism, however, prevents one from selecting only those $M$ basis functions that matter for the data at hand. This would avoid spending budget on basis functions where there is no signal, but it requires a criterion for ranking the candidates. We propose three such criteria derived from an information-theoretic view of the basis-function selection problem. Each criterion matches a different state of knowledge at selection time: a no-data state, a no-prior state, and an in-between state. We then study the performance of truncation versus selection strategies on six UCI regression benchmarks across three basis families: Hilbert-space Gaussian processes (HSGP), variational Fourier features (VFF), and variational inducing spherical harmonics (VISH). We observe that the no-data criterion is a safe default, matching or improving on truncation for HSGP, VFF and VISH, with substantial gains for VISH and improvements over a recently developed selection heuristic for that basis family. The data-aware no-prior and in-between criteria provide substantial gains over truncation specifically for HSGP, which is the most broadly used of the three families in practice.

In sparse random matrices, spectral outliers (eigenvalues and singular values located away from the bulk) emerge due to degree fluctuations: high degrees inflate the operator norm, while low column degrees reduce the least singular value. As proved by Feige and Ofek (2005) and Le, Levina, and Vershynin (2017), degree regularization enforces concentration at the expected norm scale. However, precise bounds incorporating the cutoffs remain unexplored and challenging since regularization introduces dependencies among entries. For the first time in the literature, we provide variance- and cutoff-dependent bounds for extreme singular values and eigenvalues of regularized inhomogeneous random matrices. In the absence of regularization, our lower bound for the least singular value matches the same leading constant obtained by Brailovskaya and van Handel (2024). Moreover, our error term vanishes under the milder condition $d/\log N\to\infty$, as opposed to their stronger requirement $d/(\log N)^4\to\infty$. A key ingredient is to extend spectral radius bounds for nonbacktracking matrices to the dependent setting. We build on approaches for independent cases established by Benaych-Georges, Bordenave, and Knowles (2020), as well as Dumitriu and Zhu (2024), and carefully handle edges traversed only once. Our proof framework separates deterministic spectral comparisons from probabilistic estimates: once Loewner inequalities and columnwise variance controls are established, the remaining probabilistic analysis boils down to verifying the graph moment conditions formulated in this paper. We hope this framework can be extended to handle general random matrices with more complex dependencies.

Let X1,...,Xn be independent exponential random variables with arbitrary positive, not necessarily equal, rates, and let Y1,\ldots,Yn be independent and identically distributed exponential random variables. We prove that, for every order statistic, the homogeneous order statistic is smaller than its heterogeneous counterpart in the convex transform order. Our results show that, relative to the homogeneous benchmark, heterogeneity stretches upper quantiles more strongly than lower quantiles, while the homogeneous order statistic ages faster in the convex-transform sense. We further extend the comparison to a proportional-hazards family that includes Weibull, Lomax, and Burr XII distributions. We give explicit conditions under which the comparison is preserved and show that the shape restriction is sharp within the Weibull family. These results characterize how component heterogeneity changes the distributional shape of order statistics beyond effects on location or scale.

Why do political parties polarise when most voters remain in the centre? Moving beyond conventional point-based models, we address this classic puzzle by developing a spatial voting framework that conceptualises political parties as continuous ideological distributions. Operating under endogenous turnout, these parties adapt their positioning to maximise a combination of total votes and parliamentary representation. We demonstrate that this theoretical formulation reveals a distinct phase transition in multiparty competition, moving from centripetal convergence to centrifugal separation. These dual objectives do not necessarily reward the same directional party movement as voters only cast ballots for ideologically acceptable parties. The numerical results map the precise boundaries of this transition: a higher party count destabilises the centrist equilibrium and expands the zone of divergence, whereas pronounced internal party dispersion and a stronger emphasis on representation contract it. Furthermore, electoral institutions critically govern this divergence: when voter preferences vary across geographic districts, the same national electorate can trigger entirely different phase transitions depending on how votes are converted into seats. Together, our results show how party differentiation emerges endogenously from the dynamics of electoral competition itself. Whether small initial differences grow or decay is not merely a reflection of voter preferences, but an outcome dictated by party goals, competitor count, internal party breadth, and vote-to-seat conversion rules.

In archaeological chronology building, stratigraphic constraints arising from the excavation context are commonly incorporated into Bayesian models through order-constrained priors. However, these models suffer from well-documented artifacts: estimated durations are inflated, and the bias worsens as more samples are dated. First, we prove that this phenomenon is intrinsic to the prior: when all true ages coincide, the posterior duration converges to the prior support's width. We then propose an unconstrained posterior that does not include stratigraphic knowledge and project it onto the space induced by the stratigraphic order. The projected posterior is the Wasserstein-closest constrained distribution; it can be computed via isotonic regression on a directed acyclic graph and corrects duration inflation unlike constrained priors. Simulation studies and two archaeological case studies demonstrate substantial bias reduction and precision gains.

Streaming principal component analysis (PCA) seeks to recover a leading spectral subspace in a single pass over a data stream. We give a new analysis of the ubiquitous Oja's algorithm [Oja82] for the most general, gap-free variant of this problem, where no eigengap assumptions are made on the underlying mean matrix, complemented by a nearly-matching lower bound. Prior works achieving near-optimal rates for streaming PCA either required gap assumptions [JJK+16, HNWW21], or were limited to rank-one updates [AZL17, Lia23]. Our proof only uses a second moment bound on the individual stochastic updates, bypassing the almost sure bounds needed by prior near-optimal analyses, and the analogous offline matrix Bernstein bound. We also extend our result to a Rayleigh quotient notion of approximate PCA, addressing an open question of [JJK+16]. As our main application, we give gap-free differentially private PCA guarantees for sub-Gaussian data, settling Conjecture 1.1 of [Bro26] up to logarithmic factors.

Estimating the Hüsler-Reiss precision matrix is a fundamental problem for statistical inference in multivariate extremes. In high-dimensional settings, the number of unknown parameters grows quadratically with the dimension, making regularisation indispensable. Existing approaches regularise the estimation problem by exploiting sparsity. In this paper, we consider an alternative structural assumption, namely that the Hüsler-Reiss precision matrix is block-structured. To estimate such a matrix, we introduce a new regularisation framework based on a convex fusion penalty. By encouraging rows and columns to merge, this approach provides a parsimonious representation of the precision matrix, allowing for the simultaneous estimation of its coefficients and the underlying partition of the variables. The resulting convex optimisation problem is solved by an efficient algorithm combining gradient-based updates with progressive fusion steps. We establish non-asymptotic concentration bounds for the empirical weights entering the penalty and prove consistency of both block recovery and precision matrix estimation under suitable regularity conditions. Numerical experiments demonstrate that our methodology accurately recovers the latent block structure while accurately estimating the Hüsler-Reiss precision matrix across various configurations, illustrating the practical benefits of fusion-based regularisation for multivariate extremes. These benefits are also demonstrated by applying the proposed method to foreign exchange data.

We propose a novel and flexible nonlinear approach for dimensionality reduction of large-scale multiview network data and derive its theory. The (linear) predictor incorporates observed covariates (edge-specific, layer-specific, and global) and Gaussian latent factors. Inference is conducted via the graph Laplace approximated maximum likelihood estimator. Letting $K$ denote the number of network layers and $n_V$ the number of nodes, we derive asymptotic theory under two regimes: (i) $K \to \infty$ with fixed $n_V$, and (ii) double asymptotics $K, n_V \to \infty$, establishing consistency and asymptotic normality for local and global parameters, with distinct convergence rates. In an application to the gravity model for commodity trades, we use a zero-adjusted Gamma distribution with latent factors and observable covariates (e.g.\ distance, tariffs, common language) to capture excess zeros, skewness, and unobserved heterogeneity. Synthetic and real-data exercises show that our approach outperforms the routinely applied Poisson pseudo-maximum likelihood estimator with fixed effects. We complement our theoretical and empirical contributions with open-source R/C++ routines and a novel strategy for starting values selection.

This manuscript formalizes the most popular model validation tools used in general insurance actuarial modeling. These include graphical tools like calibration plots, actual-vs-expected plots, lift charts, Murphy diagrams, as well as classical statistical tools such as Bregman losses, deviance losses, elementary losses, Murphy's decomposition and Gini scores. Particular emphasis is placed on whether calibration and discrimination are studied under a policy-weighted or an exposure-weighted population measure. This distinction is crucial in ensuring that premium schemes are calibrated on the correct scale.

This paper reviews advances over the last decade in the field of computational structural reliability modelling. The discussion is focused on problems of time-variant/time-invariant reliability at the component level. Three broad classes of problems are considered: (a) time-invariant reliability problems involving static problems, (b) problems of time-variant reliability analysis for deterministically parametered dynamical systems driven by random excitations, and (c) time-variant reliability analysis of randomly parametered dynamical systems subjected to random excitations. Four classes of approaches are discussed: (a) analytical methods based on first/second order reliability analyses (for time-invariant reliability analysis) and level crossing based approaches for time-variant reliability problems, (b) methods based on importance sampling strategies (including the Girsanov transformation based method for dynamical systems), (c) particle and trajectory splitting based methods, and (d) methods that employ machine learning based tools (primarily involving development of surrogate models and active learning strategies) in tackling reliability problems. The focus of the discussions is on methodological advances, and questions related to specific applications are not addressed. The review presents critical discussions on the relative merits of alternative approaches (in terms of computational efficiency, accuracy, scalability, treatment of rare events, ability to handle geometric complexities linked to failure surface, and suitability for high-dimensional problems) and identifies several directions for future research.

Principal component analysis (PCA) is a fundamental tool to reduce the dimensionality of the data in many applications. PCA finds a few signal directions that contain most of the variability of the data by computing the eigenvectors of the sample covariance matrix. In this work, we focus on the spiked covariance model, in which the data vectors are defined by a few orthogonal signals plus an isotropic Gaussian noise, and our goal is to estimate one or more of the leading signals. Our main theoretical finding is that the subspace spanned by several leading eigenvectors of the sample covariance matrix contains significant information about the desired signals long before the individual eigenvectors converge to the population principal components. To prove this, we derive a posteriori bounds for the angle between the subspace spanned by the desired population signals and the subspace obtained from the sample using perturbation theory for singular vectors. This leads to a new algorithm, SuperPCA (SUbsPace subsamplER PCA), which capitalizes on an approximate eigenspace of the sample covariance matrix to find the leading signals far more efficiently and accurately than classical PCA in the high-dimensional, multi-signal setting. SuperPCA exploits only a small number of subsampled coordinates of the data, which can lead to tremendous savings in data acquisition cost, especially when the signals are approximately sparse. For the same number of measurements, SuperPCA can offer a factor $10$ improvement in accuracy compared to the classical PCA method.

Longitudinal data, often in the form of panel data, are commonly encountered in the pharmaceutical industry. However it is not always possible to fit an unstructured covariance matrix when using the Mixed Model for Repeated Measures (MMRM). Here we argue that high-order autoregressive models provide a useful framework when more general correlation structures, such as the unstructured and/or Toeplitz, are unfeasible. The autoregressive model of order one is often suggested but, in our experience, higher orders are seldom considered. We present the case for using autoregressive models, with potentially high orders, as a flexible framework to model correlations in panel data, when more general structures are hard to identify.

Measuring disease burden is an important part of both public health research and health economics. Some measures of disease burden, such as hospitalisations, are made up of recurrent events. Assessing what events are related to a particular disease is however non-trivial. We extend the notion of relative survival to recurrent events by developing a novel non-parametric estimator. The estimator combines data from some cohort with aggregated population level data to estimate the number of excess recurrent events. Using empirical process theory, we show that the estimator converges weakly to a mean zero Gaussian process under mild regularity conditions, and provide a consistent estimator for the covariance function. We also evaluate the finite sample properties of the estimator through simulations and provide a practical example using data from Swedish patients with rectal cancer.

We study preference elicitation under the Bradley-Terry-Luce (BTL) model where the true partworth vector is unknown and has to be estimated as a parameter with elicited preference information. The set of selected pairwise queries is non-uniform, deterministic, and arbitrary over a collection of alternatives, provided that it satisfies a joint identifiability condition. We focus on understanding when the canonical maximum likelihood estimator (MLE) is finite and admits sharp error bounds without explicit compactness constraints on the feasible set or external regularizers. To this end, we derive minimax lower bounds under the standard bounded dynamic range condition, and find that the same Fisher-information geometry in the classic Cramér-Rao lower bounds underpins the finite-sample difficulty of the estimation problem. By combining a non-asymptotic expansion of the likelihood score equation with a fixed-point localization argument, we identify a design-dependent sample size threshold above which the unconstrained canonical MLE exists and is unique with high probability. The same expansion yields a decomposition of the estimation error into a linear stochastic term, an explicit second-order bias, and a higher-order remainder. A refined analysis gives sufficient sample size conditions under which the canonical MLE attains the minimax rates up to logarithmic and constant factors. These results provide a unified non-asymptotic theory for parametric utility elicitation and reveal when the inference is determined by response data alone rather than by external regularization. Preliminary numerical results are consistent with the theoretical findings.

Cylindrical data pair an angle with a linear measurement. Clustering that ignores the periodicity of the angle breaks up groups lying across its origin. We formulate the K-means algorithm for a generic distance on the cylinder and instantiate it with the chordal distance of the ambient space and the geodesic distance along the surface, so that the two versions differ in the metric alone. Centroids are exact Frechet means computed with known tools: the mean direction for the chord, the circular Frechet mean for the arc. Seeded by K-means++, the algorithm converges in finitely many iterations at a cost comparable to classical K-means. Over an extensive empirical analysis on simulated data, Euclidean K-means is never meaningfully better and breaks down when clusters cross the origin; the two versions agree for concentrated angles, the chord incidentally prevailing only when the angle is uninformative; a model-based cylindrical mixture is the most variable method but describes elongated, correlated clusters better. In hue-value color image quantization, the cylindrical versions performs better on object detection based on color segmentation; on wind direction and nitrogen oxides data they recover two pollution regimes that Euclidean K-means cuts apart.

Forecasters often score the same units per date against one standardized realized outcome. We show that every standardized forecast splits exactly into a component aligned with this common target and a component uncorrelated with it. Three consequences follow: forecast-error correlation largely mirrors forecast correlation and is therefore a poor measure of diversity; an equally weighted combination beats a no-information forecast only when average alignment is large relative to the combination's dispersion; and the gain from adding a forecaster separates into genuine improvement and mere dilution, which equal-weight admission can mistakenly reward. We develop a cautious selection rule, study it in simulations, and apply it to language-model forecasts of US equity rankings and mechanical signals ranking exchange-traded funds. Selection removes most dilution losses, but no combination beats the no-information forecast.

Causal tabular foundation models amortize effect estimation across synthetic mechanisms, but latent-effect supervision rewards posterior shrinkage rather than encoding the repeated-sample response needed in a fixed deployment population. We introduce fluctuation-supervised pretraining (FSP): each synthetic table is labeled by its average treatment effect plus its efficient influence-function fluctuation; deployment remains a frozen forward pass. Along the path $T_{\lambda,P}=\theta(P)+\lambda P_n\psi_P$, we prove an endpoint transition: every fixed $\lambda<1$ retains label ambiguity of order $(1-\lambda)^2/n$, whereas full fluctuation makes the Gaussian label observable and reduces optimal finite-stratum causal label-prediction risk to order $n^{-2}$. A finite-pretraining bound combines label, network, episode-sampling, and optimization errors; its sampling defect controls fixed-mechanism bias, mean squared error, variance, Gaussian approximation, and, with variance-head accuracy, studentized coverage. Complementary lower bounds separate local $n^{-1}$ ATE risk from the $\log N/M$ excess risk of generic finite-dictionary episode learning. Experiments trace the learned sampling response. Across 24 nonlinear continuous-covariate cells at trained context lengths, continuous-row FSP lowers checkpoint-mean macro RMSE by 7.0% versus S-learner and wins all 12 weak-overlap cells; validation-selected Summary FSP deploys $11.6\times$ faster per table in our warm one-thread benchmark. Under effect shift, matched Raw FSP lowers mean-checkpoint RMSE by 54.2% and teacher defect by 99.0% versus latent-effect supervision, and RMSE by 10.2% versus the released CausalPFN-S checkpoint. Known-effect semisynthesis tests coverage; two randomized-study evaluations show that lower RMSE can coexist with residual attenuation.

Explanations are central to causal reasoning, and cognitive science has long established that the human drive to explain is itself a mechanism for learning about causality. Despite this, learning from those abductive signals is largely ignored in artificial intelligence. While explainable AI (XAI) increasingly draws on causal models to generate explanations, the converse direction about what explanations can do for causality remains largely unexplored. To fill this gap, we propose xWhyL, a formal framework connecting causality and XAI by learning causal models from explanations. We develop a mathematical theory that translates explanations into a learning signal complementary to observational data, and demonstrate how it enables overcoming the limits of observational causal discovery. As explanations can be derived from incorrect beliefs and clash with data, a tension we call the Causal Tug-of-War, we prove conditions under which our framework rejects misspecified explanations rather than absorbing them. Our practical instantiation, Causal Interactive Learning (CIL), shows how expert explanations can efficiently support causal discovery and distinguish correct from incorrect explanations.

Causal inference guides operational and managerial decisions but remains challenging in high-dimensional panel or tensor settings, where decisions may depend on the joint conditional distribution of missing control outcomes. We develop Counterfactual Tucker Diffusion (\CFTDiff), which integrates the treatment mask and latent Tucker structure into conditional diffusion to recover this distribution given observed control outcomes through efficient nonlinear score learning in a low-dimensional core. The masked Tucker score preserves dependence across tensor modes while reducing the dimension of nonlinear score learning from the product of mode dimensions to the much smaller product of Tucker ranks. We establish high-probability error bounds for conditional score estimation that depend on the Tucker ranks, largest mode dimension, and the factor-strength-adjusted number of missing outcomes, and show how these bounds translate into recovery guaranties for the conditional distribution of the missing control outcomes. Across missing rates, simulations show more accurate point recovery than common causal panel and matrix/tensor completion methods; comparisons with nested diffusion specifications further demonstrate the gains from masked conditioning and Tucker dimension reduction. In Norway's iFlex experiment, \CFTDiff recovers missing outcomes more accurately than competing methods; when applied to causal analysis, its estimated conditional distributions yield counterfactual prediction intervals and target-attainment probabilities, allowing pricing interventions to be evaluated by demand-reduction magnitude and reliability.