Papers

A mixture of spatial factor analyzers (MSFA) is introduced to address the challenges of clustering high-dimensional spatial data. The proposed framework utilizes the given coordinate system to implement a flexible, spline-based spatial decay covariance structure. In this approach, the number of covariance parameters increases linearly with the number of spatial locations, rather than increasing quadratically. To model non-spatial dependence, matrix variate factor analyzers are employed for further dimensionality reduction. Parameter estimation is conducted via a variant of the expectation-maximization algorithm combined with a generalized least squares estimator. The proposed models are explored in the context of tensor-variate data analysis, where simulation studies and three applications—Raman spectroscopy, SpecTex hyperspectral textures, and the Salinas hyperspectral scene—illustrate their ability to recover and distinguish spatial dependence patterns. The R implementation is available on CRAN: https://CRAN.R-project.org/package=SpatMix.

Bubeck, Li and Nagaraj conjectured that, for generic data, any two-layer neural network with $m$ neurons that fits $n$ noisy labels must have Lipschitz constant at least of order $\sqrt{n/m}$, with no restriction on the size of the weights. Bubeck and Sellke proved a universal version of this law for Lipschitz-parameterized classes, but under a polynomial bound on the parameters; at depth three that boundedness hypothesis is genuinely necessary. The two-layer unbounded-weight case requires a different argument. We prove the conjectured law, up to one logarithmic factor, for every continuous piecewise-linear activation, in particular for ReLU networks. For data drawn uniformly from $\mathbb{S}^{d-1}$, $d\ge3$, or from $N(0,I_d/d)$, labels in $[-1,1]$ with noise level $\sigma^2>0$, and any width-$m$ two-layer network with arbitrary real weights, biases and affine skip connection, fitting the data $\varepsilon$ below the noise floor forces $\mathrm{Lip}(f)\ge c\,\varepsilon\sqrt{n/(\bar m\log(C\bar m nd/\varepsilon))}$, $\bar m=(K-1)m+1$, with high probability. A realized-kink-count version holds on the same event: every realized two-layer piecewise-linear function with $k(f)\le n$ distinct kink hyperplanes obeys the bound with $\bar m$ replaced by $k(f)+1$, irrespective of how many redundant hidden units parameterize it. The proof replaces parameter-space covering, impossible for unbounded weights, by a function-space covering. The central deterministic ingredient is a rigidity lemma: on $B_2$, and on $\mathbb{S}^{d-1}$ for $d\ge3$, the coefficient of each canonical kink is controlled by the Lipschitz constant of the realized function, because kinks on distinct hyperplanes cannot cancel at generic points. Rigidity genuinely fails at $d=2$, and an explicit two-layer ReLU interpolant with $O(1)$ Lipschitz constant at width $2n$ matches the law at the overparameterized endpoint.

Missing values undermine statistical inference and machine learning pipelines, yet most imputation methods rely on heuristics or restrictive parametric assumptions that ignore the joint data distribution. We recast imputation under missing completely at random (MCAR) as density estimation from masked observations: estimate a distribution whose observed marginals exactly match those in the data. Leveraging positive semi definite (PSD) kernel densities we obtain a convex empirical risk problem with closed form marginals, solvable by a Newton interior point method. The resulting PSD Impute model yields both single and multiple imputations from the same fitted density, enjoys statistical consistency with fast adaptive excess risk beating the curse of dimensionality for very regular probabilities. Preliminary experiments on one synthetic and eleven real world datasets already indicate competitive distributional accuracy compared with popular imputation baselines, suggesting strong practical promise.

While accuracy, robustness, and calibration are all essential for reliable neural networks, they are often studied separately; developing models that satisfy all three simultaneously remains a central challenge. Lipschitz-constrained models guarantee robustness by design, yet the manual selection of the Lipschitz constraint L governs the resulting accuracy-robustness trade-off, and their calibration properties remain largely underexplored. In this work, we highlight a theoretical and empirical link between the enforced Lipschitz constraint and Temperature Scaling, a state-of-the-art calibration method. Specifically, we find that for a given training scheme, there exists a non-trivial value L* that yields an out-of-the-box calibrated network, and that calibration acts as a principled criterion to select a well-defined operating point on the accuracy-robustness Pareto front. Leveraging these insights, we introduce Lipschitz Scaling Training (LiST), a novel training paradigm that iteratively adjusts the global Lipschitz constant to reach this operating point. Through a margin parameter in the training loss, LiST further enables the construction of a fully calibrated Pareto front, allowing users to navigate the accuracy-robustness trade-off while remaining calibrated throughout. At convergence, LiST also enables the reintegration of calibration data into training, improving sample efficiency without sacrificing calibration. We validate LiST on CIFAR-10/100 and Tiny-ImageNet, demonstrating competitive accuracy and robustness against constrained and unconstrained baselines, while remaining calibrated out of the box. Code is available at GitHub.

Phase II dose-ranging trials often report the largest observed dose-control effect while inspecting accumulating data repeatedly. This creates two coupled distortions: selection optimism from choosing the empirical winner, known as the winner's curse, and Type I error inflation from multiplicity across doses and interim looks. We develop an anytime-valid procedure for testing whether the best true dose effect exceeds a clinically meaningful margin. The mathematical starting point is a recent selection-premium identity for the running maximum: for dose-control scores, the expected gain from re-selecting the current leader becomes a predictable selection charge. Subtracting this charge gives a residual with nonpositive drift under the composite null; applying a one-sided mixture-exponential construction then yields an e-process and hence an anytime-valid global test. The resulting rule has a transparent ledger form: raw best effect minus selection charge minus monitoring margin, and a “GO” decision is made only when the remaining evidence still exceeds the clinical margin. We give plug-in implementations for Gaussian and binary outcomes, prove finite-sample Type I control and anytime lower confidence bounds for the best dose effect, and illustrate the method through a worked example and simulations.

We study minimum-norm interpolation (MNI) in overparameterized linear regression with isotropic Gaussian covariates, in settings where the MNI has no closed-form formula. Whereas most prior work relied on Gaussian comparison tools such as the convex Gaussian min–max theorem (CGMT), our approach uses tools from high-dimensional geometry and probability. First, when the norm is in isotropic position, we obtain an “offset” bound that controls the amount by which the MNI shrinks the ground truth. Second, we show that the “intrinsic” variance of the $\ell_1$-MNI is at most $O(\tfrac{1}{n\log(d/n)^2})$, using a variant of Talagrand's $L_1$–$L_2$ inequality due to Cordero-Erausquin and Ledoux [2012], together with a classical result of Gluskin [1988]. We recover the sharp mean-squared error (MSE) bound for the $\ell_1$-MNI obtained by Wang et al. [2022], using the work of Fleury [2012] on the symmetric Gaussian polytope, which is defined via \[ P_{n,d} := \mathrm{conv}\{\pm X_i\}_{i=1}^{d} \text{ where } X_i \overset{\mathrm{i.i.d.}}{\sim} N(0,\mathrm{I}_{n \times n}), \] rather than CGMT. Our methods also imply improvements on previous results in high-dimensional geometry that may be of independent interest. First, we show that with overwhelming probability, the ratio between the isotropic constant of $P_{n,d}$ and that of the Euclidean ball in $\mathbb{R}^n$ is at most $1+O((\log(d/n))^{-2})$, improving a result of Klartag and Kozma [2009]. We also establish a refined weighted thin-shell estimate on $P_{n,d}$, and provide an elementary proof of the main theorem of Fleury [2012].

Many machine learning models are defined for inputs of different sizes, such as point clouds containing different numbers of points, sequences of tokens of different lengths, and graphs on different numbers of nodes. Such models are trained on finitely many examples of necessarily limited sizes. How well do these models generalize from inputs of small size to larger inputs of size not seen during training? Furthermore, evaluating such models on large inputs is often expensive. How can we sketch large inputs to obtain smaller ones on which the model takes similar values? At the heart of both questions is the need to compare inputs of different sizes and to approximate large inputs by small ones. We present a unified approach to address these questions by using random sampling maps to compare inputs of different sizes. The sampling maps we consider are generalizations of sampling with replacement, random binning, and species sampling. We characterize the application domains in which each type of sampling is appropriate in terms of the symmetries and relations between problem instances of different sizes in the domain. Our framework yields explicit generalization and sketching rates for function classes continuous with respect to a chosen notion of sampling, encompassing large families of functions defined on sequences, graphs, and tensors of different sizes. Specific examples include moment polynomials on measures, homomorphism densities and numbers of graphs, permutation-invariant transformers, and graph neural networks.

Single-cell molecular characterization remains a bottleneck in scalable biological analysis because of labeling requirements, limited multiplexing, and reagents that perturb physiology. Raman spectroscopy addresses these limits by providing chemically specific, label-free vibrational fingerprints, but long acquisition times and specialized instruments restrict high-throughput use. Here, we overcome this barrier by showing that spectral fingerprints can be reconstructed from brightfield microscopy using generative modeling. We introduce Pic2Spec, a framework that learns a shared latent biochemical representation linking image morphology to vibrational spectral structure, enabling virtual Raman spectroscopy without hardware. We validate Pic2Spec across mammalian and bacterial cells, generating high-fidelity spectra that reproduce measured Raman fingerprints with 98% cosine similarity and Pearson correlations of ~95%, while preserving biochemical peaks and population distributions. Beyond spectral similarity, Pic2Spec provides molecular-level resolution in bacterial systems: generated spectra discriminate mutation-driven transgenic states and predict GFP expression with accuracy approaching true Raman measurements, outperforming conventional image analysis by 20%. These findings establish Pic2Spec as a first demonstration of chemically informative virtual molecular fingerprinting from brightfield images, complementing slow, hardware-intensive spectroscopy with computational inference. By redefining microscopy as an inference-enabled molecular profiling platform, Pic2Spec democratizes label-free biochemical phenotyping and overcomes the hardware and time constraints that have confined spectroscopy to specialized laboratories. This enables high-throughput molecular analysis for clinical diagnostics, screening, and monitoring at the scale and accessibility of standard microscopy.

Spatial heterogeneity in insurance risk modelling is often represented using coarse areal structures, which can obscure fine-scale patterns critical for accurate risk assessment. This study introduces a point-referenced Bayesian framework to model claim occurrence and severity at the policyholder level, avoiding reliance on predefined geographic aggregation. Drawing on a large French insurance portfolio combined with high-resolution environmental variables, rainfall records, and institutional hazard maps, we compare a benchmark GLM with several discrete Bayesian specifications, including independent random effects, intrinsic conditional autoregressive (iCAR) and Besag-York-Mollie (BYM) models, and a continuously indexed Gaussian random field constructed using the stochastic partial differential equation (SPDE) approach. Inference is performed using Integrated Nested Laplace Approximation (INLA), enabling efficient estimation of latent spatial fields and non-linear covariate effects. Our results show that accounting for spatial dependence substantially improves occurrence modelling, while gains in severity prediction are more limited. The SPDE formulation further outperforms areal models by capturing sub-municipal risk gradients and reducing artefacts induced by arbitrary geographic partitioning. By conditioning on detailed building-level attributes, we isolate the contribution of latent spatial effects, refine the interpretation of observed covariates, and improve the allocation of risk premiums across the portfolio. In addition to enhanced predictive performance, the framework provides coherent uncertainty quantification and supports tail-risk assessment. To our knowledge, this is the first application of point-referenced SPDE models to flood insurance, offering a scalable statistical alternative for pricing and managing risks with strong spatial structure.

This paper proposes a general and unified framework for testing the equality of a broad class of parameters, defined as a smooth function of expectations of symmetric kernels, across multiple independent populations. We consider two test statistics, a Wald-type statistic and an ANOVA-type statistic. The asymptotic distribution of the first one is derived under a fixed-dimension regime, whereas the second one is studied under both fixed and increasing-dimension regimes, where the parameter dimension diverges with the sample size. Based on these limiting distributions, we construct test procedures enabling asymptotically exact inference without parametric assumptions. Additionally, an alternative null distribution estimator based on a weighted bootstrap approximation is studied, which is applicable to the ANOVA-type statistic under a fixed-dimension regime. The finite-sample performance and computational efficiency of the proposed procedures are evaluated through an extensive simulation study and a real dataset application.

Harmful pesticide effects exceeding specific protection goals (SPG) may go undetected in underpowered experimental designs. Regulatory honeybee field studies have consistently failed to reach the statistical power required under European Food Safety Authority (EFSA) guidance, which may have caused approval of high-risk substances. Therefore, EFSA advised a shift from testing the null hypothesis of 'no effect' to equivalence testing. Under this approach, a pesticide is classified as 'low risk' if the null hypothesis that its effect exceeds the SPG can be rejected. For honeybees, the recommended SPG is a colony size reduction below 10%. Critics have argued that this framework requires excessive site replication to demonstrate pesticide safety and proposed an alternative equivalence test defining treatment effects relative to the lower bound of the 90%-control-group confidence interval. Using simulations mimicking a regulatory honeybee field study, we show that although the two equivalence tests share the same trade-off between false 'low-risk' and false 'high-risk' classifications, only EFSA's original recommendation reliably identifies pesticides with effects > SPG at alpha = 0.2. Our results show that increasing site replication beyond the current practice is unavoidable for a reliable regulatory assessment. However, for pesticides with effect sizes of 5% or less, site requirements remain lower than those implied by the power requirement of the former EFSA guidance. Moreover, covariate adjustment through a model term or balanced colony allocation using anticlustering randomisation can reduce site requirements without losing power and thus save costs. Finally, we provide guidance and R functions for anticlustering randomisation and equivalence testing for pesticide risk assessment.

We propose a novel approach to the minimax estimation of high-order functionals from the perspective of quantum computing. Specifically, for any real number $\alpha \gg 1$, we present two estimators, one for the classical functional $\mathrm{F}_\alpha(P) = \sum_{i=1}^S p_i^\alpha$ of a discrete distribution $P$ and the other for the quantum functional $\mathrm{F}_\alpha(\rho) = \operatorname{tr}(\rho^\alpha)$ of a mixed state $\rho$. These functionals have close connections with the Rényi entropy and the Tsallis entropy. We show that both estimators achieve the minimax optimal $L_2$ rate $\alpha \mathsf{n}^{-1}$ in the range $\alpha \lesssim \mathsf{n} \lesssim \alpha^{3-o(1)}$, where the support size $S$ of $P$ or the dimension of $\rho$ can be much larger than the number of samples $\mathsf{n}$. As a result, both estimators achieve the optimal sample complexity $\mathsf{n} \asymp \alpha$, improving upon the prior best upper bounds $O(\alpha^2)$ established by Jiao, Venkat, Han, and Weissman (IEEE Trans. Inf. Theory 2017) for classical functionals and Chen and Wang (COLT 2025) for quantum functionals. Our estimators are constructed under a unified framework using quantum primitives and run in linear time on a quantum computer. This work reveals an unexpected path from quantum computing to statistics, suggesting a conceptually new methodology for functional estimation. It adds to the growing list of quantum proofs for classical theorems.

Training a model with noisy gradient descent can be idealized as overdamped Langevin dynamics, and a natural safety question is to bound the probability $\nu_t(\mathcal{A}_H) = \mathbb{P}(Q_t \in \mathcal{A}_H)$ that the trajectory lies in a designated failure region $\mathcal{A}_H$. We study this for a smooth, strongly convex loss in $d$ dimensions, with $\mathcal{A}_H$ separated from the minimizer by an energy gap. At the end of training, the equilibrium mass $\pi(\mathcal{A}_H)$ is exponentially small in $d$, with a complementary energy-barrier rate when the noise is small. Along the trajectory, a shape-free bound $\nu_t(\mathcal{A}_H) \le \pi(\mathcal{A}_H)(1 + \sqrt{\chi_0^2/\pi(\mathcal{A}_H)}\,e^{-mt})$ shows the in-set probability relaxes to (twice) the static value after a burn-in of order $d$, using only the global spectral gap $m$. A worked Ornstein-Uhlenbeck example shows this burn-in is necessary: an angular slice of the equilibrium shell can transiently swell by a factor exponential in $d$, though its equilibrium mass is tiny. To rule this out we introduce a local relaxation rate, defined through the spectral measure of the region's centered indicator rather than a Dirichlet-form Rayleigh quotient. For geometrically isolated regions this rate exceeds the global one, shrinking the burn-in, and with a maximum-principle ceiling it caps the trajectory probability uniformly in time. Strong convexity sets how fast training relaxes, but the shape of the unsafe set decides whether the trajectory bulges through it on the way to equilibrium.

Artificial intelligence (AI) is a double-edged sword: while it has achieved remarkable success across a wide range of domains, its deployment also calls for effective oversight and regulation, for which the detection of AI-related content and artifacts is perhaps the most direct and cost-effective approach. To this end, we propose a unified detection framework based on Mahalanobis distance scores (MDS), applicable to several important settings, including the detection of large language model (LLM) generated text, hallucination, watermark, and adversarial examples. A key component of the proposed method is to accurately characterize the positive class–such as human-generated text, factual statements, unwatermarked text, or non-adversarial samples–which requires an efficient and robust estimator of the covariance matrix of deep representations of positive samples before computing the MDS. Since the positive samples typically consist of multiple classes, and these classes may exhibit both homogeneity and heterogeneity, we develop joint estimation methods for both the casewise and cellwise minimum covariance determinant (MCD) estimators. We provide efficient optimization algorithms for both estimators and prove their convergence. We provide a reasonable definition of the breakdown point for the joint estimators and prove their corresponding high breakdown point properties. Empirical evaluations confirm the effectiveness of the proposed detection framework.

We address the problem of efficiently sampling multimodal probability distributions, where standard Markov Chain Monte Carlo methods often suffer from poor mixing and mode trapping. To mitigate these issues, we propose Gradient-free Riemannian Langevin Sampler (GRiLS), a novel proposal that improves exploration without requiring gradient evaluations of the target density. Our approach introduces a Riemannian metric which reshapes the local geometry in order to facilitate transitions across modes. The resulting gradient-free MCMC algorithm is particularly suitable for complex, computationally expensive targets where derivatives are unavailable or impractical. The GRiLS proposal requires knowing the mean and covariance of the target density, which we estimate using an ensemble of interacting particles. Empirical results on multimodal benchmarks demonstrate that GRiLS achieves improved mixing compared to existing gradient-based and gradient-free MCMC approaches.

We consider the estimation of linear functionals of the mixing distribution in a nonparametric empirical Bayes framework. Our main interest is in situations in which the mixing distribution is only partially identifiable, as may arise in complex sampling situations with nonresponse. We argue that estimating the functional by applying it to the semiparametric maximum likelihood estimator of the mixing distribution is an efficient tool, even when the maximum likelihood estimator is not unique.

Self-supervised learning matches supervised accuracy from a fraction of the labels, but the labeled-sample efficiency behind this has lacked a theoretical explanation. We provide one. Data augmentation induces a similarity graph on the unlabeled data, so downstream learning on that graph is graph-Laplacian-regularized learning. We prove a fast transductive rate, $O(1/n_L)$ in the number of labels, in place of the supervised $O(1/\sqrt{n_L})$, by carrying the leave-one-out stability apparatus of Johnson and Zhang (JMLR 2007) over to the augmentation graph, and without the unrealistic assumptions of limit-based analyses (exact kernel, generalizing features). The bound makes augmentation quality explicit: the expected error is at most $C/n_L + R_{\mathrm{DA}}(y)$, where the data-augmentation alignment error $R_{\mathrm{DA}}(y)$ is proportional to the graph-cut mass of augmentations that cross a label boundary, so good augmentations let few labels suffice. The analysis uses a streamlined loss that drops the projector, negative-sample, and orthogonality overhead of standard objectives yet still recovers the top-$K$ ideal features in the infinite-data limit, the augmentation-kernel eigenspace studied by Zhai et al. The bound gives a mechanistic account of the accuracy-versus-label-count curve through augmentation quality, verified in a controlled model where the constants are known.

Background: Outcome reporting bias (ORB) occurs when study outcomes are selectively reported based on their results. ORB potentially undermines the credibility and validity of meta-analyses and contributes to research waste by distorting overall treatment effects. ORB can be viewed as a missing data problem in which unreported study outcomes introduce bias. Despite the serious implications ORB poses, it remains an underrecognized issue, with only a few adjustment methods available. Methods: We propose an approach that addresses unreported study outcomes in meta-analyses through multiple imputation for univariate and multivariate meta-analysis. To assess the impact of ORB in meta-analyses, we apply our proposed methodology to real clinical data affected by ORB, and conduct a simulation study to evaluate the method's performance under a range of scenarios. Results: The proposed method provides bias-adjusted estimates under assumed selective non-reporting mechanisms. In the application to clinical data, ORB-adjusted estimates were systematically shifted towards less extreme treatment effects compared with naive analyses, highlighting the potential magnitude of ORB in practice. The simulation study shows that the extent of adjustment depends on the assumed selection mechanism and the degree of heterogeneity, with stronger selection leading to larger adjustment. Conclusions: Imputing unreported study outcomes provides a promising approach to address ORB in meta-analyses. The multivariate approach extends ORB adjustment to jointly model correlated outcomes, allowing borrowing of strength across outcomes. Overall, we propose a practical and flexible approach for evaluating the sensitivity of univariate and multivariate meta-analytic conclusions to ORB.

We study the recovery of sparse functions from finite, noisy, and indirect observations in the framework of statistical inverse learning. The unknown is modeled as an element of $\ell^1$, and observations are generated through a possibly nonlinear forward operator $A:\ell^1\to H$, where $H$ is a vector-valued reproducing kernel Hilbert space. We propose an $\ell^1$-regularized empirical risk minimizer and develop a theoretical analysis of its statistical properties. Under mild assumptions, we establish almost-sure consistency and derive non-asymptotic high-probability convergence rates in both the prediction and $\ell^1$ reconstruction norms. The rates depend on the source smoothness parameter $r$, characterized by a variational source condition, and the effective dimension exponent $b$, describing the polynomial spectral decay of the covariance operator. We further prove matching minimax lower bounds, showing that the obtained convergence rates are optimal. To relate the theory to practical sparsity models, we consider finitely smoothing operators of the form $A=G\circ S$, where $S$ is a synthesis operator, and show that approximation-space assumptions imply the required variational source conditions. In particular, we prove that membership in the approximation space $k_t$ is equivalent to polynomial decay of the best $n$-term approximation error. Finally, we verify the assumptions for two representative inverse problems: reaction coefficient identification in elliptic PDEs and sparse computed tomography. For filtered Radon transforms, we derive explicit effective-dimension asymptotics, yielding concrete convergence rates for standard image models and sparsifying systems.

We propose an infinite hidden Markov model for cylindrical time series with von Mises-Gamma emissions. Posterior inference is performed using a beam sampler combining conjugate updates and approximate sampling schemes. Simulation studies and two real data applications demonstrate the effectiveness of the proposed methodology.

We prove that, in the realizable PAC setting, the sample complexity of exact-trace learning for full autoregressive Chain-of-Thought traces is upper bounded by the standard multiclass rate of the local next-token class, where this rate is governed by the Daniely–Shalev-Shwartz dimension. Under exact-trace loss, one wrong action makes the whole trace incorrect; nevertheless, for every stopping rule $\mathtt{halt}$ and every pointwise $\mathtt{halt}$-halting local class $\mathrm{H}$, $n_{\mathrm{PAC}}^{\varepsilon,\delta}(\operatorname{Roll}_{\mathtt{halt}}(\mathrm{H}))=O((\operatorname{DSdim}(\mathrm{H})+\log(1/\delta))/\varepsilon)$, with no dependence on rollout length. The dependence on $\operatorname{DSdim}(\mathrm{H})$ is worst-case optimal, since one-step stopping recovers ordinary multiclass learning of $\mathrm{H}$. The proof introduces parity dimension, a rollout-stable refinement of DS dimension based on even pseudo-cubes. It controls one-inclusion density via a low-coordinate spanning theorem on finite restrictions and, unlike DS dimension itself, does not increase under autoregressive rollout. We also show why this detour is necessary: DS dimension can increase under rollout.

High-order interactions are increasingly recognized as a hallmark of collective dynamics in complex systems. The relationship between high-order behaviours (HOBs), observed as synergistic or redundant statistical dependencies, and high-order mechanisms (HOMs), related to the structural or dynamical rules of the data-generating process, remains difficult to establish from data. We introduce a predictability-based framework to disentangle these two levels of description in complex network systems. Structural synergy is defined as the excess predictive power gained when two sources are considered jointly beyond the best additive description and is estimated through polynomial regression by comparing a model with interaction terms against an additive model. Simulations show that dependencies among sources reflecting HOBs can arise even in the absence of HOMs, while synergy due to a non-additive mechanism may remain hidden when the observed synergy-redundancy balance is dominated by redundancy and become detectable only through structural synergy. Applications to climate and source-reconstructed cortical EEG dynamics reveal significant non-additive predictive components despite predominantly redundancy-dominated HOBs. These findings emphasize that HOBs and HOMs can dissociate: a system may display redundancy-dominated HOBs while still containing a significant synergistic mechanism. The proposed framework supports a mechanistically informed interpretation of complex-system dynamics and may help to identify when mechanism-based models are needed to predict the response of a system to perturbations or interventions, while also recognizing that structural synergy should be interpreted as evidence of non-additive predictive structure rather than as a direct identification of the underlying generative mechanism.

Subnational life-table death counts are highly correlated across time and space and differ by gender. While these associations are helpful in improving forecasts through joint modeling, less attention has been paid to identifying and understanding mortality disparities in gender and regional gaps. We propose a forecasting framework to model and forecast female life-table death counts at the national or subnational level, and to model and forecast the associated gender gap. For either females or males, we could forecast national life-table death counts and the regional gap relative to national data. By combining gender and regional gaps, we also explore the double gap by prioritizing national data and female data. Life-table death counts are unique due to their non-negativity and summability constraints. To address the constraints, we apply a one-to-one transformation, termed cumulative distribution function transformation, to obtain one- to 15-step-ahead forecasts. Using Japanese age-specific life-table death counts between ages 0 and 110+ from 1947 to 2023, we evaluate and compare point and interval forecast accuracy across gender, region, and double gaps. By focusing on these gaps, we can deepen our understanding of the possible factors driving gender and regional mortality variations.

Accurate PM2.5 prediction depends on information available when a prediction is issued. Using 420,768 hourly observations from 12 monitoring stations in Beijing between March 2013 and February 2017, we compare same-hour estimation (h = 0) with strict one-, six-, and 24- hour-ahead forecasting. Stations are reindexed to complete hourly calendars before lagged and rolling features are constructed. Chronological train, calibration, and test blocks and expandingwindow preprocessing prevent target-time leakage. Persistence and seasonal-naive baselines are compared with Ridge, Elastic Net, and a five-seed XGBoost ensemble. Ridge attains RMSE 20.15 at one hour, XGBoost 54.13 at six hours, and Elastic Net 80.83 at 24 hours, improving on persistence by 5.7%, 11.8%, and 17.6%, respectively. Same-hour XGBoost reaches RMSE 14.03, but feature ablations attribute most of this performance to recent PM2.5 history and synchronous PM10. Moving-block intervals, repeated seeds, and leave-one-station-out evaluation quantify uncertainty and spatial transfer. The results support sensor cross-checking at h = 0, routine updates at h = 1, alert preparation at h = 6, and cautious daily planning at h = 24.

We consider the problem of minimizing the $k$-th order partial derivative $f=\partial_j^k g$ of an unknown function $g$ along a fixed coordinate direction $j$, based on noisy queries of $g$. Assuming that $g$ has Hölder regularity ${\beta+k}$ for some $\beta\ge 2$, that $f$ is strongly convex on a compact convex set $\Theta\subset\mathbb{R}^d$ and that $g$ and $f$ satisfy mild boundedness and Lipschitz regularity conditions on $\Theta$, we propose a kernel-based estimator of $\nabla f$ and analyze the projected stochastic gradient algorithm driven by this estimator. We obtain a non-asymptotic upper bound on the optimization error of the order $d^{(2\beta+k-1)/(\beta+k)}\,N^{-(\beta-1)/(\beta+k)}$, where $N$ is the total number of queries. We also establish a minimax lower bound of the order $N^{-(\beta-1)/(\beta+k)}$ showing that this rate is optimal in $N$ over all sequential algorithms.

Evidence remains limited on how missing-data strategies affect the stability of clinical prediction models across different predictor-outcome relationships and degrees of missingness. We conducted a simulation study using a fully observed real-world cardiac cohort of 8,245 patients, equally divided into development and external validation cohorts. Missing data were induced under a missing-at-random mechanism across 18 scenarios varying by variable type, predictor-outcome relationship, and missingness proportion. Five strategies were compared: complete case analysis, multiple imputation by chained equations with fully conditional specification, multiple imputation using predictive mean matching, missForest, and k-nearest neighbours. Logistic regression models were developed using backward stepwise elimination. Outcomes included optimism-corrected AUC, calibration slope, mean absolute prediction error, external validation performance, and computation time. When missingness involved isolated linear or categorical variables at 30%-60%, all methods maintained discrimination comparable to the complete-data model, with median AUCs of about 0.75. When missingness involved isolated non-linear variables or more complex patterns, predictive performance and calibration worsened as missingness increased, especially at 90%. In complex scenarios, multiple imputation showed greater prediction instability and overfitting, while missForest performed well internally but overfitted externally. k-nearest neighbours showed the most consistent performance, with stable predictions, better external validation results, and the shortest computation time. The optimal strategy may depend on the characteristics of variables with missing data. In sufficiently large development samples, k-nearest neighbours may provide a computationally efficient alternative.

We propose the Randomized BH procedure, a generalized framework that introduces randomization into the Benjamini-Hochberg (BH) procedure for false discovery rate (FDR) control. We show that the classical BH procedure as well as two recently developed families of null-distribution-free methods, competition tests and conformal tests, share a common theoretical foundation and arise as special cases of this framework. The Randomized BH relaxes the classical conditions required for FDR control, frees conformal tests from the strict exchangeability assumption while allowing calibration sets to be constructed via pre-screening, and enables competition tests to be analyzed through the well-developed theory of p-values. We illustrate the framework through two applications. The first is a deeper analysis of compound p-values through the lens of randomization. The second is a general approach for integrating distributed multiple testing problems while maintaining global FDR control. Numerical simulations confirm the validity of our theoretical results.

Diffusion models represent a leading paradigm for graph generation, with notable impact in domains such as molecular design. Yet, scaling these models to large graphs remains an open problem. We approach this question in the dense-graph setting through the lens of graphons, the size-agnostic limit objects of dense graph sequences, to study how structural graph statistics behave across node-size scales. This perspective leads to DiPhon, a diffusion framework for size-scalable graph generation. Specifically, we formulate a continuous diffusion process on the graphon space via a Jacobi stochastic differential equation (SDE), and propose DiPhon, a discretized graph-level process that mimics these dynamics on finite graphs. We further derive the corresponding reverse-time process, which requires access to the marginal score. For the Jacobi process, this score interestingly admits a tractable form, which we estimate from data via graph denoising and plug into the reverse process to generate graph samples. We prove that DiPhon matches exactly the first moment of the marginal distributions induced by the continuous graphon process, and approximates the second moment up to a closed-form discrepancy. Thus, DiPhon inherits key size-agnostic statistical properties of the graphon dynamics, providing a principled route toward scalable graph generation. Empirically, we demonstrate this scalability by training on small graphs and generating progressively larger graphs at inference time, without retraining, while preserving their core topological properties.

The Moran statistic, and its accompanying local statistics, are one of the most extensively used exploratory spatial data analysis tools for assessing global and local spatial autocorrelation. The paired visualizations for these statistics, the Moran Scatterplot and LISA map, are likewise central to spatial analysis. Together, these statistics and visualizations are used to identify spatial clusters, regions of a map where observations are similar to one another, or spatial outliers, observations that differ sharply from their surroundings. However, the use of Moran statistics to detect spatial outliers is complicated by their high sensitivity to *distributional* outliers: observations that are extreme relative to the overall data distribution, regardless of their spatial context. Indeed, a single distributional outlier can (I) distort local statistics across the entire map and (II) bias the global estimate of spatial association. Recent work has begun to address (I) and (II) separately using plug-in robust estimators for location, scale, and spatial correlation. In this paper, we offer the first systematic evaluation of robust LISA and global spatial association measures, using variety of plug-in robust estimators, a trimmed least squares (TLS) estimator, and a Theil-Sen-style estimator. We also outline a visualization strategy to create Robust Moran Scatterplots/LISA maps for each. Out of all considered approaches, we find that the Theil-Sen Moran estimator is a better default for exploratory spatial data analysis and visualization, while robust plug-in estimators also offer acceptable performance in large datasets.

Trans-model Markov chain Monte Carlo (MCMC) algorithms are widely used in Bayesian inference, and are particularly important in Bayesian phylogenetics where phylogenetic trees represent different statistical models. While the algorithm allows great flexibility, its mixing efficiency can vary hugely, and is poorly understood. Here we use mathematical analysis and simulation to explore the mixing efficiency of trans-model MCMC proposals, including the model-proposal probabilities and the proposal kernel for model parameters. Our analysis confirms the intuition that one should preferentially propose models with high posterior probabilities, and propose parameter values from the posterior as much as possible. Our results provide guidelines for constructing efficient trans-model MCMC algorithms. The principles are applied to MCMC algorithms in phylogenetic reconstruction using two real datasets for primates and mammals.