Papers

Conditional distribution equality testing is a core problem with applications in causal inference, model diagnostics, and transfer learning. Existing methods often lose power when the conditioning variable is high-dimensional. We propose a new test that leverages generative models to compare conditional distributions using both observed data and synthetic samples. A sample-splitting and studentization scheme is introduced to control estimation error. Under the null, the test statistic is asymptotically normal, and we establish local power against contiguous alternatives. Our theory under the null and fixed alternative requires only consistency of the learned conditional distribution, allowing for a broad class of generative models. Empirical results demonstrate accurate size control and strong power in high-dimensional settings compared to representative kernel-based conditional and joint two-sample tests.

Standard competition-based multiple testing methods – such as the knockoff filter and the target-decoy procedure – are typically formulated as head-to-head competitions between a test unit and its assigned negative control unit. Motivated by real-world applications, which leverage any available negative controls, we study this competition-based framework in more intricate settings where competitions vary in size and structure. We propose new methods based on compound p-values – p-values that are superuniform only on aggregate – that enable powerful multiple testing with broader competition structures, via the Bonferroni and Benjamini-Hochberg procedures. The former provides exact familywise error rate control, whereas the latter may incur a false discovery rate inflation factor, for which we derive upper bounds across a range of settings. We showcase the performance of our general competition-based testing methodology through extensive simulations and real-data applications, including analyses of virtual spatial transcriptomics and CRISPR screen datasets that exhibit precisely the complex competition structures that motivate our methods.

Regression-based survey weighting connects model-based population estimation with design-based inference through equivalent weights. For nonlinear models, equivalent weights can be defined locally through the sensitivity of a poststratified population estimator. However, conventional regression models may inadequately capture nonlinear relationships and complex interactions across sparse poststratification cells. We develop a survey-weighting framework based on hierarchical Gaussian process (HGP) regression, combining structured group-level effects with kernel-based smoothing to share information across related cells. We derive Gaussian posterior prediction formulas that exploit the structured covariance representation to reduce computational costs. We also establish a uniform convergence result for the posterior predictive mean under stated regularity conditions and construct local equivalent weights by matching the first-order sensitivity of a plug-in poststratified estimand to that of a weighted average of fitted response means, with both derivatives taken with respect to the latent mean-scale vector. Simulations show that HGP achieves lower population-weighted prediction error under nonparametric structure, whereas logistic regression performs well under correct parametric specification. Additional analyses examine sample size, frame coverage, subgroup prediction, and population-mean agreement. An application to the Future of Families and Child Wellbeing Study illustrates the method. Overall, this work extends local equivalent-weight construction to nonparametric modeling and provides a computational framework and conditional theoretical results for survey weighting in structured, sparse populations.

We study testing for peer effects in linear panel models with individual and time fixed effects when the interaction network is unspecified. We propose a quadratic-form test based on the Gram matrix of residualized peer covariates. A key difficulty is that fixed-effects removal and nuisance estimation destroy the centering of a conventional leave-diagonal-out quadratic form under conditional heteroskedasticity. We introduce a leverage adjustment that restores exact conditional centering and combine it with a residual-based variance estimator. The resulting studentized statistic is asymptotically standard normal as both panel dimensions diverge, without imposing any restriction on their relative growth rates. Under additional conditions, we establish consistency against alternatives. The null theory is further extended to covariates that are weakly dependent over time; the same studentization remains valid without an additional long-run variance estimator. Monte Carlo experiments evaluate the finite-sample size and power of the proposed test across different settings and demonstrate its favorable performance relative to the existing benchmark test. Applications to state tax competition and growth spillovers among OECD countries illustrate the practical relevance in testing for peer effects without specifying the interaction network.

Many multivariate time series arise on networks where relationships between variables depend on physical or more general distances. Standard vector autoregressive models struggle in such settings due to their quadratic parameter growth, while generalised network autoregressive (GNAR) models offer a parsimonious alternative by explicitly incorporating network structures. Existing GNAR formulations assume fixed weights, typically uniform or reciprocal distances, which may not reflect the true influence of distance in real systems. We introduce a novel model with general inverse distance weights (GNARID) that extends GNAR by introducing a tunable exponent parameter that controls the strength of distance penalisation. We establish theoretical results on the consistency and asymptotic normality of new GNARID estimators, including an extension to the setting where the number of nodes grows. We analyse the effect of misspecifying the exponent on forecasting error, which helps assess how much GNAR models could be improved by our tunable weights. Simulation studies confirm the theoretical properties and show measurable benefits from correct weight specification. Using the Brumadinho ground displacement dataset, we show that proper tuning of the weight exponent improves predictive performance, highlighting the practical value of flexible distance-based weighting in network time series models.

Fractional Brownian motion (fBM) is a model of anomalous diffusion parameterised by a Hurst exponent $H \in (0,1)$ and a generalised diffusion coefficient $D_H > 0$. We jointly estimate $H$ and $D_H$ using maximum likelihood (MLE) from density measurements with Poisson counting noise and single-particle trajectory with Gaussian localisation error. We show that one density measurement in time cannot infer both $H$ and $D_H$ and that several time points are needed, which matters most for inference precision. For trajectories, we show that a neural network trained with a squared error loss approximates the posterior mean of $H$, which can be computed exactly at negligible cost. Paired with the profile-likelihood estimate of $D_H$, the posterior mean has the lowest root mean squared error in both $H$ and $D_H$, except near $H=0$ and $H=1$ where MLE estimates $H$ better. In vesicle trajectories from two cell types, $H$ and $D_H$ form a continuum rather than distinct modes of motion. The posterior probability that each window within a track is subdiffusive or superdiffusive shows that $26\%$ of MRC5 fibroblast and $21\%$ of HeLa vesicle trajectories contain both subdiffusive and superdiffusive windows, against $2\%$ and $7\%$ expected from estimation uncertainty alone. Intracellular cargoes therefore switch between subdiffusion and superdiffusion over time rather than forming populations with fixed $H$.

Principal Component Analysis (PCA) is a standard method for dimension reduction and the exploratory analysis of multivariate data, but a complete applied workflow often requires more than numerical eigenvectors and static plots. This article documents a set of routines in the FSDA toolbox for MATLAB that implements an integrated workflow for interactive and robust PCA graphics. The software contribution consists of coordinating classical SVD-based PCA, dynamic biplot scaling controlled by row and column parameters, score-distance and orthogonal-distance diagnostics, hard-trimming robust PCA through the Minimum Covariance Determinant and the Forward Search, linked brushing, geographical displays, and three-dimensional projection diagnostics. The implementation returns reproducible numerical output while also exposing graphical controls for changing biplot geometry, confidence regions, displayed labels, robust trimming levels, and linked selections. The methodology is illustrated on quality-of-life indicators for Italian provinces, where the software reveals variable groups, atypical provinces, robust and non-robust PCA structures, and geographical patterns. The emphasis is on statistical software: no new PCA estimator is proposed, but established PCA and robust-PCA methodology is combined in a documented, extensible, and script-reproducible MATLAB environment.

This paper develops a principal component analysis framework for probability measures associated with non-negative real-valued random variables with structural zeros. Specifically, leveraging the associated Lorenz Curve (LC) and a suitable perspective transformation of the plane, a bijective embedding is defined to treat such measures as functional data in a Hilbert space. Once the random LC is embedded in this space, standard Functional Principal Component Analysis (FPCA) becomes directly applicable, and the diffeomorphic nature of the embedding allows the main convergence results and their interpretations to be transferred back to the original space of LC. The a.s. and $L^2(\mathbb{P})$ consistency of the estimated eigenfunctions, eigenvalues, and modes of variation is subsequently investigated, both in the embedded space and, by virtue of the embedding's diffeomorphism property, in the original LC space. Crucially, this consistency analysis accounts for the double stochasticity of the estimation problem: one level of randomness is tied to the underlying LC process itself, while the other stems from the empirical estimation of the LC from sample data. As a consequence, consistency of the Karhunen–Loève reconstructions is also established, guaranteeing faithful approximation of individual curves in both spaces. Finally, the proposed FPCA methodology is illustrated on a real-world dataset of global export flows and compared with a naive FPCA approach.

Background: Generalized Odds (Gen Odds) provide summaries for ordinal outcomes without requiring the proportional odds assumption. However, conventional Gen Odds compares independently drawn treatment and control outcomes and targets a population-level contrast rather than an individual-level contrast. We developed an individual-level causal formulation of Gen Odds and examined its identification. Methods: We defined the Causal Gen Odds using the joint distribution of potential outcomes and identified it under conditional independence of the potential outcomes given baseline covariates. We developed a model-based estimator using nonproportional cumulative logistic regression and introduced an odds-ratio sensitivity parameter to assess departures from conditional independence. We evaluated the estimators through simulations and illustrated them using data from the ACCORD Study. Results: In simulations, the conditional-independence estimator was approximately unbiased with 95.5% confidence interval coverage under conditional independence, but bias and undercoverage increased with residual dependence. When the sensitivity-analysis estimator incorporated the data-generating dependence parameter, bias remained below 0.01 and coverage ranged from 93.7% to 96.1% across examined settings. In the ACCORD Study, the Causal Gen Odds for right-foot pre-tibial edema was 0.909 (95% CI: 0.880 to 0.940) under conditional independence and remained stable across the examined sensitivity range. Conclusions: The Causal Gen Odds provides an individual-level estimand that complements the conventional population-level Gen Odds. This framework enables estimation of this contrast while explicitly assessing its robustness to assumptions about dependence between potential outcomes.

The unadjusted Langevin algorithm (ULA) and its stochastic-gradient variant, SGLD, are workhorses of modern Bayesian machine learning, but can be fragile when faced with a target distribution possessing a non-Lipschitz drift. In such settings, the expected squared jump distance can grow rapidly, leading to numerical instability or the need for prohibitively small step-sizes. We propose the saturated unadjusted Langevin algorithm (SULA) that pairs a saturated gradient drift with a gradient-informed state-dependent diffusion, while preserving first-order consistency with the Langevin diffusion. The drift and diffusion terms are chosen to guarantee that the one-step squared update is independent of the gradient. Consequently, SULA is robust to the choice of step-size and preserves estimation and prediction quality. We provide sufficient conditions under which SULA with a step-size $h$ converges weakly to the Langevin diffusion at the standard $O(h)$ rate, the SULA chain is geometrically ergodic, and that the bias of its invariant measure relative to the target is $O(h)$. Further, we develop SG-SULA, a stochastic-gradient version of SULA that rescales the stochastic gradient using a Gaussian CDF correction under Gaussianity of the gradient noise. Numerical experiments on various examples demonstrate SULA's robustness to tuning and SG-SULA's improved predictive performance.

Individualized survival predictions are often summarized by scalar risk scores or event probabilities at selected time horizons, although clinically relevant risk may also differ in its temporal organization. We consider finite-dimensional representations of individual hazard trajectories through interpretable hazard functionals, such as peak timing, cumulative hazard over prespecified intervals, and persistence of elevated risk. A central issue is that the richness of such representations is constrained by the hazard geometry supported by the underlying survival model. We formalize this through model-induced tangent spaces and functional ranks, derive an exact decomposition of model-supported variation into directions retained and lost by a chosen representation, and introduce a model-relative notion of local completeness. We further distinguish intrinsic model capacity from the variation actually activated by covariates. These results motivate a structured extension of AFT time scaling through a positive-affine deformation on the log-time scale, yielding a power-warp mechanism with an additional temporal deformation direction. Combined with a log-logistic reference hazard, this gives a parsimonious likelihood-based working specification for which peak timing and peak height provide complete coordinates over the unimodal regime. We also establish plug-in recovery and asymptotic inference for hazard functionals. Simulations and a lung-cancer application illustrate how model-induced restrictions propagate to hazard-functional representations.

Multivariate Hawkes models are widely used to measure endogenous interactions in event data, but latent common signals can generate apparent excitation even without structural interaction. We characterize the resulting linear projection and construct stable, nonnegative Hawkes processes that share the mean intensity and Bartlett spectrum of common-signal models despite having different excitation mechanisms. This establishes a concrete form of second-order non-identification. We derive integrated third-order cumulants for Hawkes processes with Cox immigration and establish generic local identification of normalized structural parameters in two- and three-dimensional systems. An exact-equivalence Monte Carlo experiment demonstrates how third-order information can distinguish mechanisms with identical second-order structure, while revealing substantial finite-sample limitations under weak signals. For a restricted parametric family, we develop Whittle estimation from binned counts, establish consistency and asymptotic normality under appropriate identification conditions, and construct feasible sandwich standard errors accounting for estimated mean intensities. Simulations assess structural recovery and inference accuracy. An application to cryptocurrency transactions illustrates the sensitivity of inferred excitation to model specification and the limitations of structural interpretation based on spectral fit alone.

Functional brain aging is often characterized by differences between people of different ages, although aging unfolds within individuals. This distinction is especially important for functional connectivity, where the direction of covariance change may evolve across adulthood. We introduce a Bayesian longitudinal covariance-operator framework that models each visit covariance as a latent positive-definite state informed by repeated resting-state fMRI runs. The approach calibrates run-level information for temporal dependence and replicate variability, separates within-person aging from between-person age differences, and uses an adaptive tangent representation to learn the low-dimensional structure of longitudinal change without relying on pooled covariance variation. The longitudinal operator varies smoothly in magnitude and orientation with age, while exact Fréchet derivatives translate the fitted trajectory to the covariance scale. This framework allows us to examine how brain covariance organization changes within individuals, whether longitudinal aging follows the population age gradient, and which functional systems contribute to their divergence. In the Dallas Lifespan Brain Study, within-person covariance aging differs substantially from the cross-sectional age gradient through early and middle adulthood, with the discrepancy diminishing later in life. Covariance reorganization is strongest around midlife, and the longitudinal signal is concentrated in a few coordinated directions, with most off-diagonal change occurring between functional systems. These findings demonstrate that population-level signatures of brain aging do not necessarily reflect individual aging trajectories and highlight the importance of modeling longitudinal covariance change directly.

In many modern applications, obtaining response variables is costly, time-consuming, or technically challenging, making it impractical to measure outcomes for all available units. In such settings, covariates may be available for a large pool of candidates, but responses can only be collected for a small, carefully selected subset, making the choice of which units to measure a consequential decision. Most existing subsampling methods are developed for massive datasets, emphasizing the trade-off between statistical accuracy and computational cost, and typically relying on the assumption of a correctly specified model. These features limit their suitability in measurement-constrained settings, where statistical efficiency is paramount and model misspecification is often unavoidable. To address these challenges, we develop a new subset selection framework for linear regression that accommodates potential model misspecification. The proposed approach targets predictive expected loss and yields an A-optimality-type criterion that depends only on covariate information, requiring no response observations at the selection stage. We refer to this method as Expected Loss Optimal Selection (ELOS). By focusing on predictive risk rather than model-based optimality, ELOS remains well-defined and robust under misspecification. Simulation studies and real data applications demonstrate that ELOS achieves strong predictive performance and compares favorably with existing methods across a range of scenarios.

We study inference for continuous linear functionals of solutions to ill-posed inverse problems, including Euclidean parameters arising in nonparametric instrumental variables, proximal causal inference, and shadow-variable models. Although pointwise root-$n$ inference is possible under standard conditions, we show that this need not translate into locally reliable inference. Under a model imposing only minimal regularity conditions, the target is Hellinger discontinuous when both inverse problems are ill-posed: arbitrarily small changes in the data-generating law can produce arbitrarily large changes in the target, ruling out informative locally uniform inference. Imposing a uniform bound on a minimum-norm solution restores Hellinger continuity under regularity conditions, but does not by itself ensure stable root-n inference. If, in addition, both inverse problems are ill-posed and the solution whose norm is uniformly bounded is nonunique at a law $P_0$, then under regularity conditions the target fails to be pathwise differentiable; consequently, asymptotically linear estimators exhibit nonvanishing root-$n$ bias under some local alternatives, causing Wald confidence intervals to undercover asymptotically. Under regularity conditions, pathwise differentiability is restored if the uniformly bounded solution is unique at $P_0$. It is also restored, regardless of uniqueness, under the stronger requirement that the corresponding solution satisfy a uniformly bounded $\beta$-source condition for some $\beta$>0. We relate these regimes to very weak, weak, and strong identification in instrumental variables models.

Fairness audits with several sensitive attributes are often performed one coordinate at a time. This can miss purely intersectional disparities: a score may be independent of each sensitive coordinate while depending on their joint profile. We study distributional demographic parity for a real-valued score using one-dimensional Wasserstein geometry. We distinguish two complementary constructions. Conditional and sequential audits characterize joint parity through their zero sets, but their positive values do not generally decompose a fixed joint Wasserstein disparity. By contrast, full-cell conditional quantile functions define a Hilbert-valued random variable whose variance is the joint disparity. Conditional expectations of this quantile process induce a monotone information game: every reveal order yields an exact martingale decomposition of the same joint disparity, while Shapley values provide an order-free attribution across sensitive attributes. Under independent sensitive attributes, the game reduces to a Hoeffding-ANOVA decomposition and Shapley splits each interaction energy equally among its members; under dependence, orthogonality is lost and explicit formulas reveal suppression and redundancy effects that marginal audits can hide. We give fixed-grid estimators and permutation calibration, and illustrate the framework on synthetic examples and ACSIncome. The resulting framework separates three questions that are often conflated: whether joint distributional parity holds, how a fixed joint disparity decomposes along an information path, and how that disparity can be attributed across sensitive coordinates without privileging an ordering.

While causal discovery methods have matured for complete data, such algorithms for handling missing values remain limited. We establish general sufficient conditions under which causal graphs remain identifiable (i.e., uniquely determined from the available data distribution) under various missingness categories, including Missing Not At Random, and provide a framework to extend correctness guarantees from causal discovery algorithms for complete data to missing data settings. Our results show that identifiability of causal graphs with no unobserved confounders is preserved when (1) missingness mechanisms do not causally interact, (2) no variable directly causes its own missingness, and (3) either faithfulness holds for missingness mechanisms, or they follow generalized linear models. Guided by these theoretical results, we introduce two asymptotically correct algorithms for Additive Noise Models (ANMs) that recover the underlying causal graph and provide consistent parameter estimates.

This paper extends the applicability of probabilistic graphical models to non Gaussian data. It establishes a characterization of the pairwise Markov property in undirected probabilistic graphical models in terms of copula densities. Building on this characterization, we derive a necessary condition for the pairwise Markov property in Student copula Markov networks. An Expectation-Conditional Maximization algorithm (ECM) is provided to facilitate the penalized maximum likelihood estimation and graph selection. The performance of the proposed algorithm is evaluated through simulation studies, which demonstrate its effectiveness in estimating model parameters and recovering the underlying graph structure.

Learning an individualized treatment rule (ITR) with many treatments requires estimating treatment-specific effects from limited observations per arm. Treatments with similar conditional responses can share information, but linear contrast models may miss nonlinear heterogeneity. We propose nonlinear supervised clustering fusion (NSCF), a neural Q-learning framework with a shared nonlinear covariate representation and low-dimensional treatment embeddings. Pairwise fusion groups treatments with similar learned contrasts, yielding an ITR and an interpretable partition. An alternating algorithm combines neural training, embedding fusion, validation-based partition selection, and unpenalized grouped refitting. NSCF lowers regret relative to a strong neural reference in nonlinear exact-sharing simulations; perturbations characterize the limits of hard sharing. In a melanoma xenograft application, NSCF improves average policy value over the neural reference. These results support nonlinear representation learning with supervised treatment fusion when treatment effects exhibit shared structure.

Reliable uncertainty quantification remains an important challenge in AI. Hierarchical Bayesian methods provide a principled framework for this by representing unknown quantities through posterior distributions. Joint inference of neural-network parameters and prior hyperparameters is challenging for conventional Markov chain Monte Carlo (MCMC) methods due to strong dependencies. In this work, we introduce the multilevel particle marginal Metropolis-Hastings (MLPMMH) method, which combines pseudo-marginal inference for the hyperparameters with coupled sequential Monte Carlo (SMC) across hierarchy of network-width truncations for conditional parameter inference. The hierarchy is induced by a trace-class neural network (TNN) prior with width-dependent decay, motivated empirically by the compressible structure observed in trained neural-network parameters. Under standard multilevel coupling assumptions, we establish that MLPMMH achieves mean-square error (MSE) $\mathcal{O}(\varepsilon^2)$ at expected computational cost $\mathcal{O}(\varepsilon^{-2})$, and verify the predicted rate on a Bayesian regression problem. To address larger-scale applications, we further introduce S-MLPMMH, a practical approximation that replaces the exact sampling components with a PCA-based hierarchy and scalable Bayesian Monte Carlo. Experiments on MNIST demonstrate competitive predictive performance and uncertainty quantification while exploiting the parallel structure of the approximation. Software has been made publicly available at https://github.com/submission-code-anon/mlpmmh-anon.

Network-linked data combine a network with covariates and responses observed for its nodes. Common regression approaches either model dependence among neighboring responses using the observed network or capture similarities associated with network structure without explicitly modeling response dependence. We propose a probabilistic network-linked data analysis framework that incorporates both features while accounting for randomness in network formation. The model uses normalized connection probabilities to weight peer responses and latent positions to capture network cohesion. A two-step likelihood procedure first estimates the network structure and then fits the response model. We establish theoretical guarantees for network recovery in sparse networks with heterogeneous node effects, and prove consistency and joint asymptotic normality of the regression estimators subject to suitable conditions. The analysis also clarifies when estimating the network preserves the classical convergence rate for spatial autoregressive models. Simulation studies demonstrate the estimation performance of the method and assess its inference accuracy. An analysis of diplomatic network data illustrates its practical application.

Design calculations let researchers gauge, before data is collected, whether a planned study is likely to detect an effect. In this paper, we develop such calculations for predictive models using ordinary least squares and ridge regression under leave-one-out cross-validation. To do so, we derive closed-form expressions relating the expected out-of-sample $R^2$ and the expected correlation between observed and predicted outcomes to four design-stage quantities: sample size, number of predictors, an assumed population $R^2$ ($\rho^2_{in}$) and a design efficiency factor $D^*$ (which summarizes how leverage is distributed across observations). Unlike the classical shrinkage formulas, which depend only on sample size and dimension, these are sensitive to how covariates are arranged – to imbalance and skew. A simulation study assesses the accuracy of the formulas against realized values under a variety of data generating mechanisms. We apply the methods to a published study predicting personality from smartphone behavior and show that the design formulas would have indicated, before any analyses were run, that least squares could not be expected to beat the sample mean while ridge could. An R package, predpower, implements the calculations.

Impact-echo testing is a nondestructive testing method used to detect internal defects in concrete structures. However, boundary reflections, structured interference, and measurement noise can produce complex Fourier amplitude spectra, making defect classification difficult. We constructed a two-dimensional frequency–amplitude persistence image (FA-PI) from the Fourier amplitude spectrum. Using a fixed parameterization of exact multi-parameter persistent homology, each retained Fourier component is represented with frequency-bin location as one coordinate and amplitude as the other, and Gaussian kernels are then applied in this two-dimensional space. This construction has a direct analytical connection to the corresponding one-dimensional Gaussian-smoothed Fourier amplitude spectrum. FA-PI was evaluated on an experimental backfill-void dataset and a public bridge-deck impact-echo dataset, and compared with raw and Gaussian-smoothed Fourier amplitude spectra using a linear classifier with a Softmax output and identical signal-level data splits. The evaluation included original signals and signals corrupted by Gaussian noise, sinusoidal interference, or mixed noise, with matching noise types and signal-to-noise ratios during training, validation, and testing. Performance was assessed using accuracy and macro-averaged precision, recall, and F1-score. Across the two datasets, FA-PI consistently outperformed the raw and Gaussian-smoothed Fourier baselines in the noisy-condition summaries. These findings indicate that representing Fourier spectral components in a two-dimensional frequency–amplitude space can improve impact-echo defect classification beyond frequency-axis smoothing alone.

We construct confidence sets for the entire Fréchet mean set by inverting studentized pairwise risk inequalities. Gaussian multiplier calibration localized around empirical near minimizers gives asymptotic coverage at least the nominal level, uniformly over distribution classes satisfying empirical-process conditions. The construction requires neither unique ness nor a manifold structure. We characterize the limiting binding directions and give conditions for exact asymptotic coverage. On the intrinsic circle, exactness holds for every fixed non-point-mass law and uniformly over a family in which two means merge at a quar tic minimum. With ideal multiplier quantiles, the construction reproduces the Euclidean Wald set under both nonsingular and singular covariance. Under regular local conditions, it is first-order equivalent to intrinsic Wald inference.

In many applications one object is chosen to represent a sample and the least typical objects are set aside, usually without any measure of uncertainty. We measure how central an object is by its reference depth, the probability that a random reference point is closer to it than to a random opponent. Using only pairwise dissimilarities, we construct simultaneous intervals for all ranks, a centre set that contains the most central observation under stated conditions, and outlier sets that contain only objects among the least central ones, each with asymptotic probability at least a prescribed level. These guarantees hold even though the ranked objects are themselves random. For symmetric distributions on the real line, the centre set has the smallest possible order of size up to a logarithmic factor, and in dimension three or more it reduces to the most central observation with probability tending to one. For typical meteorological years the centre set usually contains several candidate years, so the data do not determine the year on which building-energy simulations rest; for three collections of gene trees it contains one tree, whose features can be reported as those of the most central tree.

Recently, Morris et al. introduced TinyLoRA, a strategy for LLM fine-tuning which introduces arbitrarily few parameters, and demonstrated accuracy improvements using less than $10$ parameters rather than the ${\gtrsim 10,000}$ typically demanded by LoRA-based methods. However, that study continued to use gradient-based local search methods of the kind deployed for training multi-billion parameter LLMs in the first place. This article is predicated on the observation that, in the low dimensional setting, global optimization methods become tractable, unlocking Large language Model Bayesian Optimization (LlMBO). LlMBO is sample-efficient, gradient-free, memory-light, parallelizes across devices without synchronization, and, unlike policy-gradient methods, can optimize greedy-decoding accuracy directly. We describe numerical experiments comparing LlMBO against GRPO on fine-tuning tasks for a variety of models on a math benchmark suite which show that LlMBO can be a suitable alternative in this regime, matching or exceeding GRPO's fine tuning gains, while using a fraction of the compute and memory.

The most widely used investment risk measures - variance, standard deviation, and tracking error - treat upside deviations as equally undesirable as downside losses. In stable markets with assumed symmetric return distributions, this simplification is potentially harmless. In the climate transition age, it is not. Climate risk is two-dimensionally asymmetric, across assets and over time. Both dimensions are negatively skewed against investors. Hence, climate transition investment strategies deliberately deviate from benchmarks that reflect an insufficiently decarbonising global economy. When fiduciaries enforce tracking error budgets that penalise outperformance just as heavily as underperformance, they inadvertently constrain exactly the portfolio return opportunities that climate-aware investing seeks. We argue that replacing symmetric deviation-based risk measures with their downside-only counterparts - semi-variance, standard risk and trailing error - resolves this conflict without sacrificing any fiduciary discipline. The distinction is not merely semantic: as the gap between climate transition trajectories and market benchmarks widens over the coming decades, symmetric risk budgets will become binding constraints on transition strategies well before actual investment risk warrants it. We propose specific changes to mandate design that would allow institutional investors to pursue climate transition opportunities within genuinely risk-aware (rather than deviation-averse) frameworks.

We study a class of strategic hypothesis testing problems in which a principal aims to allocate resources among competing agents who possess data or data-generating mechanisms that may be relevant to the principal's goals. Each agent can pay a participation cost to opt into a statistical contract, in which case they generate data from a private distribution and receive an allocation based on a predetermined (and publicly known) mechanism. Each agent is either "null" or "non-null" and the principal's goal is to design the contract to limit the allocation to the former while maximizing transactions with the latter. Using e-values, we define a "continuous-closure," a set of allocations which prevents null agents from profiting in expectation and provides post-hoc risk control for general loss functions under arbitrary dependence. When the loss is the fraction of resources allocated to null agents, the region admits polynomial-time linear optimization. We also give explicit, interpretable closed-form approximations, and we analyze participation from a non-null agent's perspective, giving conditions under which entry into the contract is profitable. Finally, we extend the framework to support online arrival of the agents, in which case the principal has to make irrevocable on-the-fly decisions, and provide a method which gives time-uniform error-control in this setting.

We propose a multilevel main effects matrix factor model (MMEFM) to extract meaningful row and column effects from groups of matrix data. As a generalization of a multilevel matrix factor model, we provide rigorous definitions of MMEFM and identifications, together with iterative algorithms for the estimators of all global and local components in the main effects and common component, including global/local main effects and core ranks, with rates of convergence spelt out. An important feature of the model is that all main effects can be nonstationary and still be consistently estimated at each time point. At the same time, our estimated global common factors are alignment independent, and hence are truly global, and can be presented and given interpretations more easily. We propose a Bootstrap procedure for detecting the existence and identifying the groups of data that are driven by the global common factors. The power of the test, together with the accuracy of various estimators, are demonstrated in extensive simulation settings. Estimation performance is also compared to other state-of-the-art methods for matrix time series by simulations and a set of U.S. equity data. All numerical works are done using the R package MMEFM, which is available at CRAN.

Stochastic differential equation (SDE) models can reproduce trajectories while missing temporal dependence. We compare residual-based scores from a held-out path with scores from paths simulated by a fitted implementation, keeping fitted parameters, simulation rules and the scoring rule fixed. We give conditions for Monte Carlo calibration and bound conditional rejection error by the gap between conditional score CDFs. Residuals use zero-noise centres and covariance approximations that track noise propagation over observation intervals; report residuals also include reporting variation. Residual weights affect sensitivity and the influence of fitting error. For a quadratic energy score in a linear example with singular true covariance, whitening with nonsingular fitted covariance matrices leaves a nonvanishing score-CDF gap as nonzero coefficient estimates approach the truth. A fixed positive-definite weighting matrix makes the gap vanish. In a controlled SIR experiment using reference paths simulated at fitted parameters, adding susceptible-state endpoints to training or downweighting a low-variance residual reduces null rejection. A separate comparison using reference paths simulated at the true parameter shows reduced sensitivity to added susceptible-coordinate noise after downweighting. An exploratory weekly-report application incorporates reporting and reconstruction into reference simulations, revealing excess correlation between squared residuals in Czechia and persistently negative residuals in Finland.