Papers

Engression (Shen and Meinshausen, 2024) learns a conditional distribution by fitting a generative model $Y = f(X,\varepsilon)$ under the energy score, a strictly proper scoring rule. We provide a theoretical error analysis of engression implemented with deep neural networks. We decompose the excess risk into three components: the approximation error, the stochastic error, and the Monte Carlo error. Based on this decomposition, we establish convergence rates under the assumption that the target conditional generator admits a compositional smoothness structure.

Estimating the Shannon entropy of discrete distributions with countably infinite support is a challenging problem. In this paper, we investigate the bias of the plug-in estimator $\hat{H}_n$ for the Shannon entropy $H(\boldsymbol{p})$ under an occupancy model whose frequency sequence exhibits regular variation with tail index $\alpha\in(0,1)$. Using Poissonization and the theory of regular variation, we establish the asymptotic relation $|\mathsf{E}[\hat H_n] - H(\boldsymbol{p})| \sim n^{\alpha-1}L(n)C_\alpha$, where $L$ is a slowly varying function and $C_\alpha$ is an explicit constant depending only on $\alpha$ that admits an integral representation. Our result shows that the asymptotic behavior of the bias of the plug-in estimator under power-law frequency distributions is determined by the tail behavior of the underlying distribution.

In algorithm-assisted decision-making in high-stakes settings like healthcare, an algorithmic decision support tool provides a recommendation, but the human ultimately makes the decision. Determining whether algorithm assistance actually improves the quality of human decision-making prior to deployment is critical, and randomized experiments are one way to collect robust evidence. Historically, however, experimental designs and analyses ignore how decision-making behavior adapts with repeated algorithm exposure. In this work, we define a set of effect estimands that account for and characterize human behavior adaptation under repeated exposure and justify why these estimands are useful to target for developing a better understanding of the impact of algorithm assistance. We propose a minimax stepped double wedge design that facilitates estimating these target estimands. Finally, we compare our proposed design to two common alternative designs identified through a review of historical randomized trials of algorithm assistance. We show that these designs are less amenable to estimating the target estimands and produce biased estimates under three different forms of behavioral adaptation inspired by dynamics observed in the real world – automation bias, alert fatigue, and calibrated reliance.

Comparing graph partitions is fundamental to the analysis of network-structured data, yet existing measures for comparing graph partitions typically rely on graph-agnostic indices that treat vertices as exchangeable, ignoring the underlying graph topology that encodes essential information about community cohesion and separation. We propose a general construction of graph-informed distances that compares vertex partitions through induced edge partitions and yields valid metrics on the space of contiguous graph partitions. As special cases, we develop graph-informed versions of variation of information and the van Dongen distance together with a binary cut-based companion distance, and show that these distances satisfy a natural local graph-aware refinement criterion. Under stochastic block models, we prove that stronger topological disruptions incur asymptotically larger distances almost surely in both inter-community and intra-community split settings. These results provide a simple and principled framework to compare graph partitions while respecting the underlying graph structure.

Motivated by the study of heterogeneous returns to education in Brand & Xie 2010, which considers how the effect of completing college on earnings varies with the (unknown) probability of completing college, we analyze the problem of estimating a nonparametric regression function when certain covariates are estimated in a first step. Plug-in estimators that treat the estimated covariates as known generally suffer from first-stage estimation error. To mitigate this issue, we analyze two debiasing approaches within a framework that is agnostic to the choice of the first-stage estimation method and relies on either local-smoothing or sieve-based methods for the second-stage regression. In particular, we consider: (i) influence function-based estimators of pathwise differentiable parameters that approximate the target estimand, and (ii) a variant of plug-in estimators that directly aims to correct their bias. For each method, we upper bound the estimation error and characterize conditions under which oracle rates can be approached, highlighting the possible gains in terms of convergence rates relative to the plug-ins. Simulation studies illustrate the finite-sample behavior of the methods. We apply our methodology to data from the National Longitudinal Survey of Youth 1997 and find evidence that completing college yields the largest reductions in unemployment for individuals least likely to do so, consistent with earlier findings in the literature (Brand & Xie 2010; Brand 2023).

Survival inference often requires uncertainty quantification under censoring and limited sample sizes, while prior information may be available but difficult to incorporate without specifying a full likelihood. Likelihood-based Bayesian survival methods provide prior-informed inference but may be sensitive to distributional assumptions and computationally demanding for complex survival models. The Bayesian bootstrap provides a likelihood-free, nonparametric approach to uncertainty quantification, but by itself does not provide a general mechanism for incorporating priors on model parameters. We propose a general Bayesian method for survival models that combines the generalized Bayesian (Gibbs) updating with Bayesian bootstrap. The Bayesian bootstrap generates a distribution over a target survival estimator through Dirichlet weights, while generalized Bayesian updating incorporates parameter-specific prior information through a loss function. The resulting posterior provides a flexible alternative to likelihood-based Bayesian inference and can be applied to a broad class of survival estimators. We develop the method for the Cox proportional hazards model and obtain posterior inference for regression coefficients and hazard ratios. Simulation studies demonstrate uncertainty quantification and prior-data learning across varying sample sizes. An application to right-censored survival data illustrates its practical utility. The methodology is implemented in the open-source R package BayesBoots.

Causal mediation analysis serves as a key tool for uncovering the mediating mechanisms linking treatments to outcomes. Existing methods for mediation analysis with multiple mediators typically assume complete observations or missing-at-random and may yield biased estimation when mediator values are missing not at random (MNAR). This paper studies the identification and estimation of causal mediation effects with multiple mediators subject to MNAR missingness. We consider a broad class of MNAR mechanisms in which missingness may depend on unobserved mediators, treatment, covariates, and outcomes. Under a series of increasingly general MNAR mechanisms, we establish identified natural direct and indirect effects, effectively generalizing existing mediation analysis to handle nonignorable missing mediators. Based on the proposed identification framework, we develop estimation procedures for causal mediation effects and evaluate their finite-sample performance through simulation studies. The results demonstrate satisfactory performance across a range of missingness scenarios. An application to data from the National Health and Nutrition Examination Survey(NHANES) illustrates the practical utility of the proposed methodology for investigating mediation pathways in the presence of nonignorable missing data.

Wavelength selection is one of the important preprocessing methods in near-infrared spectroscopy to improve prediction accuracy and interpretability of spectral data. We formulate wavelength-region selection for sugar content estimation as a binary black-box optimization problem and propose a method based on Bayesian optimization. The proposed method constructs a sparse quadratic surrogate model and sequentially extracts interested wavelength regions by Thompson sampling. Minimizing an acquisition function is performed as a quadratic unconstrained binary optimization problem by simulated or quantum annealing. Experiments show that the proposed method improves the prediction accuracy of partial least squares regression and yields more consistent wavelength regions than genetic-algorithm-based selection and simulated annealing. Under one-bit local perturbations, the selected wavelength regions show minimal fluctuations in root mean square errors between observed and predicted values of a validation set. This local stability suggests that our method converges to a smoother error landscape and avoids isolated overfitted solutions. These results indicate that combinatorial Bayesian optimization is a useful framework for robust feature selection in spectroscopic prediction tasks.

The causal inference literature has traditionally focused on estimating the mean of the potential outcome, whereas evaluating how a treatment affects the entire outcome distribution can provide additional information in biomedical research. Quantile treatment effect (QTE) captures such distributional differences, particularly when outcomes are skewed. However, existing approaches for estimating QTE make distributional assumptions about the outcome and are thus sensitive to model misspecification. Motivated by an HIV study with skewed outcomes, one of which is subject to detection limits, we propose a doubly robust framework for estimating QTE based on the cumulative probability model (CPM), which is a rank-based, semiparametric linear transformation model. We develop two CPM-based estimation strategies: (1) an inverse-cumulative distribution function (CDF) approach that first estimates the marginal CDF of potential outcomes using the efficient influence function (EIF) and then obtains marginal quantiles via weighted quantile interpolation by inverting the distribution, and (2) a direct approach that solves the EIF of potential marginal quantiles. The proposed estimators are doubly robust and asymptotically normal. We further extend the framework to probability treatment effects (PTEs) and their conditional counterparts. For statistical inference, we investigate several variance estimation procedures, including EIF-based estimators, sandwich estimators, and the nonparametric bootstrap. Simulation studies illustrate that the empirical sandwich estimator and the nonparametric bootstrap provide doubly robust variance estimation with stable finite-sample performance under nuisance model misspecification. The proposed methods are evaluated through extensive Monte Carlo simulations and illustrated using an HIV data application.

Sampling representative nodes from large graphs is fundamental to graph signal processing and network analysis, yet existing methods require access to the full graph Laplacian, making them impractical at scale. We propose a simple and effective column-selective graph sampling algorithm based on a minimum inner product greedy selection rule. At each iteration, the algorithm accesses only a small random subset of Laplacian columns, requiring no eigendecomposition or global graph traversal, making it well-suited for large-scale graphs where the full Laplacian cannot be stored in memory. We analyze the algorithm under the stochastic block model and show that, when the degree distribution is balanced across nodes, the algorithm achieves sampling proportional to cluster size, and that the resulting mean estimate is controlled for band-limited graph signals in the Paley-Wiener space, with the error decaying as inter-cluster connectivity weakens. Numerical experiments on both synthetic and real-world data validate the effectiveness of the proposed method.

Unmeasured confounding is widely recognized as a limitation of observational causal inference, but its implications for data-driven causal discovery are often understated. We provide an asymptotic characterization of this limitation. When the true causal variable is unobserved but correlated proxy variables are available, increasing sample size while searching for increasingly stable, data-driven causal structures will make associations involving observable proxies become increasingly certain. This result highlights a fundamental distinction between reducing statistical uncertainty and recovering causal mechanisms: more data can improve certainty about genuine relationships within the observed variable space while providing no guarantee of identifying the underlying causal structure.

The non-Gaussian filter provides a deterministic numerical method for nonlinear and non-Gaussian state-space models, but its application has long been limited due to the computational cost of numerical integration. Advances in computing power and memory capacity have substantially reduced this limitation for low and moderate dimensional models. This paper re-examines the non-Gaussian filter and demonstrates its usefulness for maximum likelihood estimation and Bayesian inference. Numerical experiments with linear, nonlinear and radar-tracking models show that log-likelihood obtained by non-Gaussian filter is smooth and can be optimized reliably, whereas the ones obtained by particle filter are affected strongly by Monte Carlo variability even with many particles. Bayesian estimation is performed using a self-organizing state-space model, in which unknown parameters are incorporated into the state vector and estimated jointly with the latent states. These results demonstrate that the current computing technology has renewed the practical value of deterministic filtering for statistical inference in state-space models.

New venture start-ups need to “survive” through multiple stages of reaching milestone targets. We investigate the strategies for start-ups in a milestone-oriented setting. We examine a model of an entrepreneurial start-up firm, where its state is captured by a diffusion process. The entrepreneur can choose between multiple activities (or controls), which incur different cost and determine the drift and the variance of the process. Depending on whether the process reaches a fixed upper boundary or a lower one, the start-up firm succeeds or fails. Continuous-time stochastic models with multiple ($\ge 3$) controls are typically very challenging to deal with. In this work, we are able to completely solve for the optimal policy and provide an explicit characterization of its structure. In particular, the optimal policy only uses controls from a set characterized by a so-called efficient frontier curve that orders the controls by two intuitive measures: riskiness (drift-to-volatility ratio) and cost-effectiveness (drift-to-cost ratio). A unique feature of our model is that depending on the model parameters, the efficient frontier curves can be of different types, resulting in qualitatively different structures of the optimal policy. As far as we know, this is the first study that analyzes a stochastic control model which admits efficient frontier curves of different types. Our work provides start-up firms with intuitive measures to evaluate their activities and offers valuable insights on how the optimal strategies in a milestone-oriented setting change qualitatively contingent upon the specific scenario. We believe the results provide a foundational block in the study of entrepreneurial decision-making.

While randomization justifies the use of unadjusted effect estimators in randomized trials, there is growing interest in covariate adjustment to improve precision. Adjusting for baseline variables that are prognostic of the outcome can reduce estimator variance, resulting in narrower confidence intervals and increased statistical power. Recent guidance by the U.S. Food and Drug Administration supports fixed adjustment for prognostic covariates using parametric regression models. However, this guidance does not address more flexible approaches using data-adaptive or machine learning methods. We offer our perspectives on covariate adjustment to improve analytic precision. We focus on estimating the average effect for the target population in trials with minimal outcome missingness. We provide a non-technical overview of effect estimators that are unadjusted and effect estimators using fixed versus data-adaptive adjustment. We offer practical suggestions for conducting adjusted analyses that are data-adaptive, fully pre-specified, transparently and reproducibly implemented, robust to model misspecification, and guaranteed to improve precision relative to unadjusted analyses — all while preserving statistical validity and the causal effect of interest. We hope that sharing our perspectives will foster broader discussion and eventual acceptance of principled, pre-specified, data-adaptive covariate adjustment in randomized trials.

An LLM answering a question about a relational database must know which tables hold the information. Showing too few makes the question unanswerable, while showing too many wastes context. Most schema-linking methods score each table alone, though the required tables are usually connected by foreign keys and can include connectors matching the question only weakly. We develop a hierarchical Bayesian model that combines question evidence with the foreign-key graph and returns empirically near-calibrated probabilities over complete sets of tables, not a single list, so later stages can tell when a selection is uncertain. Schema graphs are small enough to evaluate every subset exactly. With a fixed context budget the graph improves the chance of retaining every required table, and the benefit grows with schema size, turning positive at about seven tables. When the model also chooses how many tables to keep, a learned prior on subset size accounts for the whole improvement and the graph adds nothing. Fit the size prior first, and use the graph when the budget is fixed. The model's predictive sets cover the true set about as often as promised, where independent scoring is overconfident.

In staggered treatment-adoption designs, later-treated units are valid controls for an earlier-treated cohort only until their own treatment begins, so the admissible donor set contracts with event time. Fixing the donor pool at the longest horizon discards temporarily eligible donors, whereas re-estimating synthetic-control weights independently at each horizon can make the counterfactual unstable as donor composition changes. We propose Risk-Set Transported Synthetic Control with Difference-in-Differences Adjustment (RT-SC-DiD). For each cohort and event-time horizon, the estimator fits weights on the currently untreated donors while shrinking them toward a transported reference that reallocates the weight of exiting donors to similar surviving donors. A DiD baseline correction removes persistent level differences. We characterize distortion from horizon-by-horizon reoptimization, derive the loading change induced by naive deletion and renormalization, and give a conditional recursive bound for error propagation under an explicitly assumed regularity condition on the transport map. We also introduce donor-support diagnostics and a donor-only placebo procedure for selecting the transport penalty. In an 80-replication pilot comparison and a separate 40-replication-per-value sensitivity analysis, intermediate transport regularization reduces average RMSE relative to independent horizon-specific estimation and strong anchoring. This evidence supports the method's bias-variance motivation but is not a proved guarantee. RT-SC-DiD is intended for settings where later-treated units provide useful short-horizon information and donor support contracts materially over time. Existing staggered synthetic-control and synthetic difference-in-differences methods do not, to our knowledge, explicitly regularize within-cohort weight sequences toward transported references as risk sets contract.

Retrieval plug-ins supply a deep forecaster with information its lookback window cannot carry. Published evaluations report consistent gains, and each credits its own mechanism. We show that the benefit belongs instead to the operating point: the relation between window length $S$ and dominant seasonal period $L$, an axis the standard protocol never varies. Stratifying the evaluation by that relation exposes the regime. At $S{=}12$, a simple control that repeats the last observed period beats the six standard backbones, in aggregate, on four of seven benchmarks by $8\%$ to $44\%$ of MSE. It beats the strongest plug-in we run on ETTm1 and matches it on ECL. It is worse by up to $25\%$ on the three datasets whose training-split spectra lack a concentrated, shared period. A controlled synthetic sweep of horizon, period, and window shows the benefit boundary tracks the period (correlation $+0.71$), not the horizon ($-0.23$). A paired control with no phase to recover nearly erases the effect, consistent with phase starvation. Zero-shot pretraining does not escape it: a foundation model trails trained backbones by $22\%$ to $50\%$ on the periodic benchmarks. Within our instrument, exact lookup matches graph diffusion: the payoff is consulting the record, not the machinery on top. Two interpretable statistics, a trend test and a staleness rate, predict the sign of the per-cell benefit at $0.76$ accuracy under leave-one-dataset-out evaluation, a suggestive margin over the $0.69$ majority rule, where a 22-feature stack manages $0.57$. We propose no new plug-in. The contribution is the regime map, the protocol that reveals it, and two statistics that screen it before deployment. Code: https://github.com/KurbanIntelligenceLab/retrieval-regime.

Modern networks can be huge with millions or even billions of nodes and edges. Thus, algorithms must be capable of scaling to such large networks in order to be practically useful. In this work, we are interested in developing an algorithm to identify core-periphery structure in massive networks. Core-periphery structure is a meso-scale feature where nodes are grouped into a densely connected core or sparsely connected periphery. To identify such structures in large networks, we propose a divide-and-conquer algorithm. The key feature of our algorithm is leveraging the edge list representation of the network, instead of the adjacency matrix, as it tends to be faster and makes a more efficient use of memory. We apply the proposed algorithm to synthetic and real-world data, notably demonstrating its performance on a real-world network with almost 14 million edges without loading the entire network into memory.

Football injury incidence is conventionally reported per 1000 hours of exposure, dividing counts by recorded playing time. Where recorded minutes track the workload under study, that division removes part of the association it should scale. We aimed to quantify the two defects contaminating a recorded-minute denominator, outcome truncation and exposure-dependent playing time, and to measure the denominator gradient across leagues and sexes. Retrospective observational cohort study of 88,573 appearances by 1208 male English Premier League players, 2017-2025, with replication from appearance records alone in 628,487 men's and 29,799 women's appearances across 15 European leagues. Changing only the offset moved the seven-day estimate from a rate ratio of 1.27 (95% CI 1.11 to 1.44) to 1.09 (0.95 to 1.25), an attenuation of 0.151 (0.140 to 0.163) that divided away 14.0% of the rate ratio; the outcome truncated recorded minutes but explained none of it (-0.095%, resampled interval -0.41 to 0.20). Independent sources confirmed attribution in 24 of 27 resolved records (88.9%, 71.9 to 96.1). The gradient exceeded the negligible threshold of 0.05 in all 15 leagues: 0.431 to 0.655 in men's competitions, 0.444 to 0.647 in women's. Restricting to starters returned it below that threshold in eight of eight men's leagues but five of seven women's, the bound reaching 0.058 in Sweden. Recorded minutes are not a neutral denominator, and the gradient should be reported before any per-minute rate. These are measurements from public records, not a tested hypothesis, and the remedy is league-specific: run the check locally.

Heavy tails weaken high-confidence control for the empirical mean. Geometric median-of-means (MOM) also lacks a threshold that moves toward mean efficiency. We propose HOMER, or Huber-of-Means for Efficient and Robust Estimation. HOMER aggregates block means through a radial Huber center. Its canonical and pseudo-Huber forms bound each block score and interpolate between median-like robustness and the empirical mean. We establish a Hilbert-space majority theorem and a MOM-order deviation bound under a finite second moment. Canonical HOMER recovers the sample mean inside its quadratic region. Pseudo-HOMER approaches the sample mean as the threshold grows. It also admits asymptotic linearity and consistent sandwich covariance estimation around the population block-Huber target. Under a finite third moment, fixed finite-dimensional projections support mean inference at the usual parametric rate. This result requires growing block sizes and counts, with block sizes increasing faster. Heavy-tailed simulations show that HOMER remains stable when a minority of block summaries is displaced. On clean Gaussian data, both versions closely approach the empirical mean's efficiency. Finite-block sandwich intervals undercovered, especially for skewed functional data. Further studies show failure when contamination affects most blocks or compromises ordinary within-block means.

Longitudinal biomedical studies increasingly collect irregularly sampled, multivariate categorical data that imperfectly reflect disease progression. This presents three key analytical challenges: highly heterogeneous disease progression across subjects; the presence of unobserved, co-evolving latent variables driving multiple measurements; and highly irregular sampling intervals both within and across patients. To address those challenges, we present the Continuous Latent Ornstein-Uhlenbeck Dynamics (CLOUD) framework for modeling complex disease trajectories from multivariate longitudinal categorical data. CLOUD links multivariate categorical observations to underlying latent functional domains and characterizes their coupled temporal evolution via the integration of a measurement component adjusted from item response theory (IRT) with a dynamic component based on multivariate Ornstein-Uhlenbeck (OU) processes. Methodologically, we introduce a time-inhomogeneous OU process that incorporated covariate-dependent components into the shifting mean function of the latent dynamics, allowing baseline biomarkers to modulate individual-level disease trajectories. We further propose a scalable parameterization of the OU drift matrix that enabled valid interaction modeling without restricting the number of latent functional domains. We establish theoretical properties of the proposed framework, including the analytical tractability of the model, the generality of the drift matrix reparameterization, and the identifiability of the entire model. Through simulation studies and an application to longitudinal amyotrophic lateral sclerosis (ALS) clinical data, we demonstrate that CLOUD provided a principled and flexible tool for characterizing subject-specific disease evolution across multiple interacting functional domains.

In this paper, we study the sample complexity of the empirical plug-in estimator for the $2$-Gromov-Wasserstein distance $D_2$ between compactly supported probability measures on Euclidean spaces. Let $\mu$ and $\nu$ be supported on compact subsets of $\mathbb{R}^{d_x}$ and $\mathbb{R}^{d_y}$, respectively, and let $\widehat\mu_n$ and $\widehat\nu_n$ be their empirical measures based on independent samples of size $n$. We prove that \[ \mathbb{E}\left|D_2^2(\widehat\mu_n,\widehat\nu_n)-D_2^2(\mu,\nu)\right| \lesssim n^{-2/((d_x\wedge d_y)\vee 4)} (\log n)^{\mathbf 1_{\{d_x\wedge d_y=4\}}}. \] This rate is sharp up to the logarithmic factor in the critical dimension. The proof is based on a geometric representation of the Euclidean distance as a squared $L^2$-distance between half-space feature maps. This yields a variational dual formulation of the Gromov-Wasserstein functional in terms of a family of classical optimal transport problems indexed by an infinite-dimensional auxiliary parameter. Although the resulting cost functions need not be semiconcave in either argument, we introduce a marginal recentering of the costs that restores the concavity structure needed for sharp metric-entropy bounds. Combining this representation with empirical-process estimates gives a rate governed by the smaller of the two ambient dimensions.

We prove an $\Omega(L/T)$ lower bound for the convergence rate of minimization in the class of functions that are convex and $L$-smooth relative to negative entropy on the standard $d$-simplex, valid for every first-order method when $d = \Omega(T^2)$. In particular, this shows that mirror descent is optimal up to a logarithmic factor in this class. This may be surprising due to the fact that accelerated methods are readily available under the assumption of smoothness in $\ell_1$-norm. While Dragomir et al. (Mathematical Programming, 2022) have already showed that acceleration might be impossible under relative smoothness, their prox-function is pathological and constructed together with the hard instance. In contrast, we show non-acceleration for a specific prox-function with particularly favorable structure. We also extend the result to the quantum setting, proving the same lower bound in the class of functions $L$-smooth relative to negative von Neumann entropy on the spectrahedron of $d \times d$ Hermitian positive-semidefinite matrices with unit trace.

Two-block data—two sets of variables measured on the same individuals, such as microbial taxa and metabolites—raise the question of how groups of covariate variables relate to groups of response variables. Co-clustering answers this for a single matrix, not for two variable blocks; existing two-block methods either cluster only one side or return signed factors rather than clusters. Starting from the multivariate linear regression $Y_1\approx M Y_2$, we give its non-negative coefficient matrix a tri-factorization $M=X_1\Theta X_2$ (a tri-NMF), so that $X_1$ softly clusters the response variables, $X_2$ the covariate variables, and $\Theta$ is a tested matrix of block correspondences. This makes the method the non-negative member of the reduced-rank regression (RRR) family, expressing RRR's low-rank class in a parts-based basis as NMF relates to PCA; the constraint can only restrict the fit, so predictive accuracy is not the aim; the co-clustering and tested correspondences are. We give multiplicative update rules, choose the two ranks by cross-validation, and develop a conditional Wald test for $\Theta$ applied after basis selection; its size is nominal with fixed bases, conservative after re-estimation, and slightly above nominal under correlated responses, while a non-zero path's magnitude stays conditional on the estimated bases. We illustrate the method—a tri-factorized non-negative RRR (NMF-RRR)—on four data sets spanning a permutation structure (Doubs, community ecology), a weak cross-structure under $p>n$ (nutrimouse, nutrigenomics), a pronounced one in a screened microbiome–metabolome study (FRANZOSA, where two microbial groups are jointly associated with each metabolite module), and a classification special case (Wine).

Gaussian Process Regression (GPR) is a robust framework for uncertainty quantification, yet its $O(n^3)$ complexity limits its scalability. Low-rank Nyström approximations can reduce this burden to $O(nm^2)$, but their accuracy depends heavily on the selection of landmark points. We propose an adaptive Nyström approach that greedily selects landmarks to minimize the trace residual of the kernel approximation error. Unlike static approximations, our method interleaves landmark expansion with hyperparameter optimization, allowing the selection process to adapt as the covariance structure is refined. Numerical experiments on five benchmark functions demonstrate that this method significantly outperforms random landmark selection in both accuracy and stability. It achieves predictive performance comparable to exact GP inference while maintaining linear scaling with respect to the sample size, providing a principled and efficient framework for large-scale computer experiments.

We develop a posterior-level discrete-to-continuum theory for Bayesian inverse problems whose finite-dimensional priors arise from local finite-difference regularization and whose likelihoods are of a general convex GLM-type form. In contrast to continuum-first approaches, the continuum prior and posterior are not assumed at the outset but are constructed as limits of the finite-dimensional probability measures used in numerical computation. First, with the total information parameter $\tau$ and the regularization strength $\kappa$ fixed, we prove weak convergence on $L^2$ of the reconstructed discrete Gaussian priors and posterior measures to well-defined continuum laws. We then consider a coupled grid-refinement and small-noise limit in which $N\to\infty$ and $\tau_N,\kappa_N\to\infty$, while $\kappa_N/\tau_N$ remains fixed. Under explicit growth conditions relating $N$ and $\tau_N$, we prove that the reconstructed discrete MAP estimates converge to the unique minimizer $u_*$ of the limiting continuum cost functional, that the reconstructed posterior measures concentrate at $u_*$, and that their centered and rescaled fluctuations converge to $\mathcal N(\vec 0,Q^{-1})$, where $Q$ is the Hessian of the continuum cost functional at $u_*$. Finally, we show that the same deterministic limit and Gaussian fluctuation law are obtained by first constructing the continuum posterior and then taking its small-noise, high-regularization limit.

Incomplete compositional data analysis faces a fundamental limitation: likelihood-based methods for compositional models generally require fully observed compositions, making it difficult to accommodate missing or censored proportions directly on the simplex. Consequently, analysts often discard partially observed compositions or transform the data into unconstrained spaces, potentially sacrificing interpretability and coherence. This paper proposes a likelihood-based method for incomplete compositional data without leaving the simplex. Specifically, we develop an Expectation-Maximisation (EM) type algorithm for fitting finite mixtures of Dirichlet distributions in the presence of missing and censored components. The proposed approach performs parameter estimation and model-based imputation simultaneously while preserving the compositional structure and interpretability of the original variables. A simulation experiment evaluates the performance of the proposed estimators and imputations under increasingly complex coarsening mechanisms. Particular attention is paid to clustering performance, and model selection outcomes. The results showed beneficial clustering performance despite observations being incomplete, and a higher probability of model selection metrics identifying the correct number of clusters compared to current alternative of case-deletion. The practical utility of the method is illustrated using two real datasets with distinct coarsened patterns. Analysis of the xenolith dataset identifies a four-component Dirichlet mixture that reveals interpretable profiles of rock types and speciation methods. Application to PM$_{2.5}$ speciation data from the Air Quality System, containing both left-censored and missing-at-random values, supports a four-component mixture model that characterises compositional parts of particulate matter across the United States.

A novel measure of dependence between a continuous random variable and a multinomial random variable is introduced. The proposed measure is based on the Hellinger distance between conditional distributions. It satisfies the desiderata for a dependence measure without making specific distributional assumptions about the continuous random variable or assuming that the discrete random variable arises from a latent continuous random variable. An estimator of the dependence measure based on data splitting and kernel density estimation is developed. The asymptotic distribution of the estimator has a simple form with a convergence rate of \sqrt{n}, making confidence intervals for the dependence measure and a test for independence straightforward and computationally convenient.

In switchback experiments with unequal cluster sizes, outcome dispersion across randomization units inflates estimator variance and limits statistical power. Standard control-using-prediction-as-covariate (CUPAC) adjustment may be suboptimal for variance reduction in this setting, because in its basic form it targets overall predictive accuracy and does not distinguish between the components of variance that vary across the randomization units of switchback experiments and those that vary across individual observations, even though these components contribute unequally to estimator variance. We propose a power-optimal variance reduction methodology via CUPAC that balances prediction of the noise between and within randomization units to achieve maximum statistical power. The methodology utilizes the framework for decomposition of the variance of the treatment-effect estimator for switchback experiments, and adapts both the outcome prediction and the analysis-time residualization to minimize the treatment-effect variance. The study first develops the theoretical framework for the power-optimal CUPAC. We then validate the theoretical framework through an extensive Monte Carlo simulation study. Finally, we discuss the practical considerations of the proposed methodology, including its potential efficiency gains and limitations.

Temporal heterogeneous graphs offer a natural abstraction for dynamic relational systems in which diverse node and relation types co-exist and evolve over time. Learning on such graphs requires jointly modeling cross-type structural heterogeneity and the temporal dynamics of interactions, yet existing methods still struggle to reconcile parameter-efficient cross-type transfer with relation-aware specialization, and typically inject time only as additive features outside the attention kernel. We propose THGFM, a web-scale temporal heterogeneous graph fusion model that addresses both limitations within a unified dual-path architecture. THGFM couples a Shared-Space Temporal Attention branch for parameter-efficient cross-type transfer with a Relational Type-Partitioned Temporal Attention branch for relation-aware specialization, and integrates them through Dual-Path Relational–Shared Fusion, instantiated with Type-Conditioned Non-Competitive Gated Sum Fusion: a adaptive mechanism that assigns independent, type-conditioned feature-wise gates to the shared and specialized branches, allowing both to be amplified or suppressed without zero-sum competition. To directly incorporate relative time into the attention score, THGFM further introduces Rotary Temporal Attention, which rotates queries and keys by half-phases of relative time before matching. THGFM consistently outperforms baseline graph transformer models on academic graphs benchmarks, delivering a $+3.25\%$ six-task mean gain, with peak relative gains of $+12.37\%$ on OAG-CS PV, $+4.87\%$ on PF-$L_2$, and $+1.18\%$ on PF-$L_1$, and $+4.24\%$, $+3.73\%$, and $+4.61\%$ on OGBN-MAG, HTAG-ArXiv, and HTAG-DBLP, respectively.