Papers

Wildfire spread poses substantial environmental and public-health risks, motivating interpretable models for short-term forecasting. We develop a statistical framework that combines cellular automata with ideas from network epidemiology to model wildfire evolution across a spatial lattice. Each grid cell is classified as available, burning, or consumed. State transitions distinguish spread from burning neighbors, intrinsic ignition, and cessation of burning, with transition rates linked to meteorological and environmental covariates. A likelihood-based estimation procedure yields transition-specific covariate effects and probabilistic forecasts of cell states. We assess forecasting performance in a simulation study and in applications to the 2018 California wildfires and the 2019-2020 Australian wildfires. We also compare the method with a published forecasting approach using the 2017 Haypress fire. The results show strong short-term discrimination in many settings, with reduced accuracy at longer forecast horizons and during abrupt fire expansion. The framework provides an interpretable basis for studying wildfire dynamics and identifies opportunities to improve ignition forecasts through richer spatial and observation models.

The problem of partial identification concerns bounding causal quantities that remain unidentifiable given the observed distribution and the causal diagram of the underlying structural causal model (SCM). While bounds have been developed for certain causal quantities with continuous outcomes, they are necessarily univariate and are unable to reflect differing levels of confidence. Building on previous work in partial identification, we propose a simple Bayesian method for deriving probabilistic bounds on a large class of causal quantities with potentially multivariate and/or continuous outcome variables.

We consider nonparametric regression with a Banach-space-valued predictor where the covariate takes values in a general separable Banach space. Our primary objective is to develop a unified framework for estimating the regression operator without relying on an inner-product structure, thereby extending classical and functional nonparametric regression beyond Euclidean and Hilbert spaces in the presence of measurement error. We establish the fundamental properties of the proposed estimator and develop its large-sample theory under conditions formulated in terms of the geometry and local concentration of the Banach-valued predictor. In particular, we study pointwise estimation and inference and discuss the construction of pointwise confidence intervals and uniform confidence bands. Our proposed framework is further extended to accommodate error-in-variables settings, where the predictor is observed with contamination. We also investigate related inverse problems arising from the construction of inverse weighting functions and consider partial contamination models. These developments provide a general framework for nonparametric regression and inference with complex infinite-dimensional predictors and identify several statistical and mathematical challenges that arise when the usual Hilbert-space structure is unavailable.

This paper introduces a novel framework for Regularized Multivariate Functional Principal Component Analysis (ReMFPCA) via Functional Singular Value Decomposition (SVD). The proposed method extends existing MFPCA approaches by incorporating a generalized functional SVD within a Hilbert space framework, enabling simultaneous regularization of both functional principal components (PCs) and their associated PC scores. A key innovation of this framework is the inclusion of a sparsity penalty on the PC scores, which enhances interpretability by filtering out irrelevant subject-specific variations. This dual-penalization strategy represents a significant advancement beyond existing covariance-based eigen decomposition methods, which penalize only the functional PCs. Two power algorithm implementations, sequential and joint, are proposed, together with a cross-validation approach based on iterative regression for optimal smoothing parameter selection. Comprehensive simulation studies and real data applications demonstrate that the proposed framework substantially improves the extraction of informative and interpretable components, offering methodological and practical benefits for analyzing multivariate functional data across diverse domains.

High dimensional spatio-temporal models can quickly run up against computational limitations. Recursive or multi-stage Bayesian algorithms are one way to address this computational challenge. Here we consider an increasingly used two-stage algorithm. The first stage models locations independently and in parallel across space. Resampling particles from this stage-one model with Metropolis-within-Gibbs methods and suitably computed acceptance ratios can impose spatial dependence across locations at a low computational cost, while targeting the same posterior distribution as the single-stage MCMC algorithm at a fraction of the overall computational cost. In this paper we show three complete examples of how two-stage MCMC algorithms can be used to efficiently fit Bayesian spatio-temporal models to large datasets. Specifically, we consider: a spatio-temporal self-exciting count model on areal data utilizing intrinsic conditional autoregressive (ICAR) prior distributions; a spatio-temporal model of point-referenced data using a stationary, isotropic distance-based covariance function; and a binary model fit to data on a spatial lattice utilizing a latent Gaussian process to capture all spatio-temporal dependence. In each example, we define the model and describe the corresponding two-stage MCMC approach. We then compare the resulting posterior samples and computational efficiency with those obtained using a standard single-stage MCMC algorithm. While varying across the three examples and numerous parameters, the two-stage approach often achieves roughly an order of magnitude higher computational efficiency, with the two-stage and single-stage algorithms producing closely agreeing posterior samples.

Many modern prediction tasks involve data from multiple heterogeneous sources, while the test distribution may differ substantially from any individual source. Although heterogeneity poses challenges, it also offers an opportunity: different sources may provide complementary information, with some regions of the feature space better represented in one source than another. We propose Multi-Source Randomly Localized Conformal Prediction (MS-RLCP), which builds on the local coverage properties of randomly localized conformal prediction (RLCP) (Hore and Barber, 2025) and extends it to multiple sources through data-adaptive source selection. Under the widely adopted assumption of a shared response distribution conditional on the features across sources and the test population, we establish finite-sample coverage bounds using an interpretable notion of envelope distribution that captures their aggregate feature-space representation. Our analysis allows the test feature distribution to be absolutely continuous with respect to the envelope, extending beyond mixtures of source distributions. Under additional regularity conditions, we also establish asymptotic test-conditional coverage. Simulations and real-world experiments demonstrate the effectiveness of MS-RLCP across varying levels of data heterogeneity.

A covariance matrix rarely distorts in a single direction: a shift can expand all variances simultaneously, shift variance along one or two principal directions, displace spectral mass without changing marginal means, or rotate the dependence structure. Default mean-shift detectors are blind to these effects, and no existing online procedure calibrates a family of spectral deviation tests with provable false alarm control across both time and tracking window selection; this paper fills that gap. We consider four spectral deviations: $D_{KL}(P_{1}\|P_{0})$, $D_{KL}(P_{0}\|P_{1})$, Jeffreys, and Bhattacharyya, each evaluated on the eigenvalues of the empirical relative covariance operator $\widehat{\Sigma}_{0}^{-1/2}\widehat{\Sigma}_{t}\widehat{\Sigma}_{0}^{-1/2}$. Critical values come from a conditional parametric bootstrap that accounts for estimation uncertainty in both the past and tracking windows, an aspect asymptotic approaches typically overlook. The procedure controls the false alarm rate family-by-family over a predefined monitoring period and a range of candidate window sizes; when a single operational window is needed, a power-based criterion selects it. We prove consistency under constant alternatives via local spectral expansions with second-order sensitivity near the null. Simulations are conservative under the null and show detection power depends largely on spectral shape rather than magnitude: $D_{KL}(P_{1}\|P_{0})$ excels under global inflation, while Jeffreys and Bhattacharyya cover a broader range of alternatives. We illustrate the approach on three financial applications: European stock indices, Fama-French sector portfolios, and large-cap technology stocks.

Many interventions act on an entire spatial exposure field. We study when their induced laws satisfy policy positivity and what its failure implies for causal evaluation through weighting and structural modelling. We show that, for post-transformation fields that are Gaussian, the only admissible translation that satisfy positivity lie on the Cameron-Martin space of the observed exposure field. Under fine-scale variation and regularity assumptions, these are the only admissible non-decreasing pointwise policies on that scale. Common rules, including binding caps and proportional reductions on the Gaussian scale can therefore violate positivity. The Cameron-Martin norm controls both weight variability and the amplification of disagreement between candidate outcome models into disagreement between policy estimates. Without positivity, bounded continuous response models can make arbitrarily similar observational predictions while their policy estimates remain separated. We also give a policy-specific criterion for structural identification conditional on a realised exposure state. Finally, when a policy violating positivity admits weights in every finite representation, their second moments diverge under nested refinement that recovers the full field. These results motivate assessing identification and sensitivity for the specified policy effect across plausible structural models and under refinement of finite representations.

In this study, we propose a novel multivariate Nakagami-$m$ distribution whose joint probability density function admits a fully closed-form expression. The proposed density is derived using a Mathai-type construction based on the partial sums of independent gamma random variables. For this distribution, the maximum likelihood estimator is not available in the closed form. To address this limitation, an asymptotically efficient closed-form estimator is developed by applying a one-step refinement to a $\sqrt{n}$-consistent initial estimator. Monte Carlo simulations demonstrate that the proposed estimator achieves a performance nearly identical to that of the numerically computed maximum likelihood estimator, while consistently outperforming the closed-form initial estimators. A real-data application further illustrates its practical utility.

Two frequent approaches for identifying structural VARs are external instruments, which carry economic content but are often weak, and non-Gaussianity of the shocks, which provides statistical identification but carries no economic meaning. We combine the two strategies in a single generalized method of moments framework. This hybrid approach point-identifies the target shocks while also identifying the non-target shocks up to sign and ordering. Under suitable rank conditions, the higher-order moments of the structural shocks anchor the identification uniformly over the instrument strength. Consequently, under local-to-zero proxy relevance, estimators of the dynamic causal effects remain consistent, and standard asymptotic inference remains valid. Moreover, the Anderson-Rubin confidence sets are substantially narrower than their instrument-only counterparts. The hybrid estimator is also more efficient than either source of identification used in isolation: at any fixed proxy relevance, even a weak instrument increases the efficiency of the estimator through its covariance with the non-Gaussian moment block. Under local deviations from proxy exogeneity, we provide asymptotic bias bounds and show that stronger non-Gaussianity of the shocks compresses the bias. Finally, the over-identified structure yields two mutually orthogonal specification tests for proxy exogeneity and the validity of the higher-order moment conditions. We derive their limiting distributions and provide a bootstrap procedure for finite-sample critical values. Monte Carlo simulations and two applications with identification of oil news-shock and a Euro-area MP shock demonstrate the potential of our framework.

High-dimensional distributions can change without altering means or covariances, while the number of affected coordinates is often unknown. We propose nonparametric procedures that address both challenges through standardized rank comparisons of marginal distributions. Sum and maximum scans target dense and sparse changes, and a Cauchy combination adapts to unknown sparsity. The procedures require no marginal moment assumptions and extend to multiple-change detection through wild binary segmentation. Under a weakly dependent Gaussian copula model, we establish asymptotic null distributions, asymptotic independence, detection consistency, and localization guarantees. Simulations demonstrate competitive performance for changes in shape and tails, including alternatives that preserve the first two moments. Applications to gene expression and sensor data illustrate the practical benefits of combining dense and sparse evidence.

When analyzing heterogeneous data with latent subpopulations, finite mixture regression (FMR) models are effective as they allow for variations in regression coefficients across mixture components. Variable selection is important at two levels: (1) at the group level by removing completely irrelevant covariates across all subpopulations and (2) at the individual level by removing irrelevant covariates within each subpopulation. However, existing variable selection methods for FMR focus solely on individual-level regression coefficient selection, providing no mechanism for group-level variable elimination. We introduce a component-weighted sparse group LASSO regularization method that performs efficient variable selection at both levels simultaneously in finite mixture regression models belonging to the exponential family. We develop a novel optimization procedure through a Majorization-Minimization algorithm. Simulation studies under Gaussian and Poisson mixture regression settings demonstrate reliable recovery of the true sparsity structure at both levels. The method is illustrated by analyzing two real datasets.

Distributed multiple testing asks $N$ sites, each holding p-values for its own hypotheses, to control the false discovery rate (FDR) of the discoveries made across the whole network while communicating only a small number of bits. Greedy interval aggregation (Pournaderi and Xiang, IEEE TSIPN, 2023) meets the communication budget but controls FDR only asymptotically, and its FDR can exceed the target at finite sample sizes; the same interval counts both select the rejection regions and calibrate the stopping rule, which biases the selected interval densities upward. We propose budgeted BONuS-GA: each site mixes synthetic uniform p-values into its data, the center ranks candidate p-value intervals across sites using the pooled counts, and each site's false discoveries are estimated from its own synthetic counts under a per-site share of the level. The procedure needs no knowledge of the null proportions, uses every p-value for both selection and inference, and controls $FDR\le\alpha$ at every sample size for any fixed assignment of hypotheses to sites, assuming only that the null p-values are independent uniforms, independent of the non-null p-values. It keeps the original $O(\sqrt m\log m)$ communication, $m$ being the total number of p-values, when $N=O(\sqrt m)$. We also give a sample-splitting alternative, cross-fit greedy aggregation, with finite-sample control under a random-site model, and we quantify the selection bias behind the original procedure's failure. In a monitoring-network simulation the budgeted procedure retains most of the power of its heuristic counterpart at moderate-to-large site sizes.

Isotonic regression and mean estimation over a simplex are among the most classical denoising problems under geometric constraints. Chatterjee and Lafferty generalized both to the recovery of a flow on a rooted tree from Gaussian noise, with isotonic regression on a path and the simplex on a star. Their results reach only special trees; on a general tree the minimax risk remained unknown for a decade. We settle the question in full. On every tree with $n$ vertices, at every budget $V$ and noise level $\sigma$, the minimax risk is of order $\min\{V^2H_T,\ \sigma^2k_0\}$, with $H_T$ the diameter and $k_0$ the crossing index of a truncated ancestor profile. A deterministic estimator, an exponentially weighted aggregate over integer states, attains it in $O(n^2\log(n+1))$ arithmetic operations. Least squares at a known budget, the natural convex program for the model, is suboptimal in the worst case by a factor of order $(\log n)^{2/5}$. One functional, the truncated ancestor profile, carries the rate, the estimator, and the least squares comparison. It serves as covering radius, as packing entropy in the positive cone, and as the range of the one-dimensional state behind the exact computation. Two classical theories become one on trees, sharp and computable throughout.

Variable importance scores from tree-based methods such as random forests favor continuous predictors over categorical ones. We present a theoretical analysis of this bias and propose a simple remedy: add a small amount of noise to each categorical predictor. The correction is demonstrated on a variety of simulated and real-world datasets and combined with integrated path stability selection to perform variable selection with false discovery control for mixed data.

In light of the US FDA announcement supporting the use of Bayesian methods in clinical trials, we present a nontechnical review of problems with Bayesian analyses and methods to address them. Our focus is on the well-known sensitivities of Bayesian results to biases embedded in their prior distributions. This vulnerability calls for special safety features, including detailed justification of informative priors based on empirical-research literature rather than on singular expert opinions. Crucially, the information contributed by priors should be realistic in light of actual background data rather than merely judged so by available experts. We illustrate our recommendations with using a textbook clinical-trial example which covered these points. We then use the example describe basic diagnostic procedures and safety recommendations for Bayesian analyses and their presentation. These include reference analyses which exclude informative priors, and justification for informative priors based on empirical-research literature rather than on philosophical arguments. We also critique some common priors that violate our realism requirements. One strategy to aid judgments about priors is to translate them into data to be added to the study data. The source of these prior data should be examined with the same critical attitude as the source of study data. We illustrate some simple summary methods for comparing prior and data information, including prior predictions of data, variance ratios, and their translation to effective sample sizes (ESS). Finally, we argue that Bayesian analyses should be treated as extensions or supplements to conventional analyses, rather than replace them, for that treatment puts prior construction and justification at the fore of the transition from frequentist to Bayesian analyses in teaching and practice.

One of the main advantages of the Bayesian approach to statistical inference is the flexibility in incorporating information from various sources, from expert opinion to historical data. Whereas the literature on Bayesian dynamic borrowing is rich, practical guidance for how to implement such methods is comparatively sparse. We thoroughly review the most common dynamic borrowing approaches for external data, including how to incorporate multiple historical data sets. We then show how these techniques can be used in practice, using publicly available software, to perform important statistical tasks such as model selection, average treatment effect estimation and prior sensitivity analysis. The aim is to make analysts more confident in their use of Bayesian methods to incorporate external data. All code is made publicly available in a companion GitHub repository (https://github.com/EzequielEBS/hdbayes-tutorial).

Developing an asymptotic theory for restricted maximum likelihood (REML) estimators in generalized linear mixed models (GLMMs) has remained an open problem for decades. In this paper, we establish the asymptotic distributions of the REML and maximum likelihood (ML) estimators for GLMMs when both the number of clusters and the cluster sizes tend to infinity. Under very mild conditions, requiring only finite-moment assumptions on the random effects rather than normality and imposing no restriction on the relative rates at which the number of clusters and the cluster sizes diverge, both estimators are proved to be asymptotically normal with an explicit block-diagonal covariance structure, and the bias-correction property of REML is established. Simulation studies support the asymptotic results, and a genetic data analysis illustrates the methodology.

The Swing filter extends a feasible line through its previous recording and restarts from a fitted value at the last accepted sample. We show that this inherited endpoint causes a singular loss of compression on rough input. For uniformly sampled Brownian motion at fixed tolerance, every fixed number of recording times converges to the first continuous feasibility boundary, with alternating boundary errors. The longest post-first segment tends to zero in probability. The segment-count ratio against a continuous approximation through true sampled knots diverges almost surely whenever that boundary precedes the observation horizon. For unit Gaussian increments, constrained least squares and no maximum lag, the first mean segment length has quadratic order in the tolerance $E$, the second has linear order, and every fixed mean from the third onward has order $E\log E$. We explain the logarithm through the conditional survival tail and prove a geometric bound on the subsequent approach to the stationary coefficient in the large-tolerance limit. The stationary mean satisfies $C_E\sim\gamma_{\mathrm{Swing}}E\log E$. An explicit one-dimensional kernel specifies the single coefficient $\gamma_{\mathrm{Swing}}\approx1.8120703$ and admits a $3/4$ Wasserstein contraction and convergent evaluation brackets. The decimal is not a certified error interval. A half-budget slope corridor, true endpoint chords and one-sample bridges repair the degeneration while preserving the original error bound, continuous output and constant working memory. Including bridge costs, the long-stream mean span per output segment is asymptotic to $21\zeta(3)E^2/(8\pi^2)$. Its explicit coefficient follows from the known Brownian anchored lifetime. Thus a change to the recording rule restores the quadratic scale.

Low-rank representations such as empirical orthogonal function (EOF) decompositions are widely used for analyzing large spatio-temporal environmental fields. However, conventional EOF identifies latent modes solely from covariance structure and does not utilize observed environmental covariates, limiting its ability to incorporate external information into low-rank representations. This study introduces Structured Covariate-Informed EOF (SCIEOF), a covariate-informed extension of EOF that bridges low-rank dimension reduction and prediction-oriented spatio-temporal modeling. SCIEOF embeds spatial and temporal covariates into the latent bases while incorporating spatio-temporal covariates through an additive component, yielding low-rank representations with latent modes informed by observed covariates. Estimation procedures are developed and evaluated through simulation studies and an application to global near-surface air temperature from the MERRA-2 reanalysis. Simulation studies demonstrate that incorporating informative covariates improves latent structure recovery and predictive accuracy, particularly when the spatial basis is appropriately specified. The advantage is more pronounced at moderate-to-large sample sizes, while methods with stronger structural assumptions remain competitive when data are limited. In the MERRA-2 application, SCIEOF achieves competitive or improved predictive performance relative to commonly used methods while providing a compact and physically interpretable low-rank representation. Overall, SCIEOF provides a flexible and computationally scalable framework for integrating structural covariate information into low-rank spatio-temporal representations, extending EOF toward predictive environmental modeling.

In this paper, we propose a novel nonparametric estimator for the correlation coefficient that is based on solving a cubic equation. This approach allows for the estimation of time-varying correlation coefficients in a nonparametric framework, providing flexibility in capturing complex relationships between variables. Furthermore, we adopt the local linear smoothing technique to correct the boundary effects, which are common in nonparametric estimation. We establish the theoretical properties of the proposed estimator, including consistency and asymptotic normality, and demonstrate its performance through simulation studies. Additionally, we apply our method to analyse dynamic functional connectivity in brain networks under six frequency bands, highlighting its potential for uncovering insights into neural interactions and cognitive processes.

This contribution develops a surface decomposition method for first-passage dynamic reliability sensitivity analysis of linear systems exposed to non-Gaussian wind excitations. The first-passage failure probability sensitivity is formulated as a system surface integral over a highly non-smooth and high-dimensional hypersurface. This complex integral is first decomposed into a collection of component surface integrals over the truncated smooth quadratic hypersurfaces. The dominant components are then identified based on the relative magnitudes of the first-order approximations of the component failure probabilities. A nested sampling algorithm is constructed to efficiently estimate the sum of these component surface integrals, in which the number of system limit-state function evaluations equals the number of outer-level samples while remaining independent of the inner-level sample size. A key advantage of the present approach is that the function evaluation results can be reused across different design parameters. Two numerical examples are explored to demonstrate the effectiveness of the proposed method. The results indicate that the number of function evaluations required is typically below 100 to achieve a target coefficient of variation of 0.1.

Staircase cluster randomized trials restrict outcome collection to a limited window around each cluster's transition from control to intervention. This can reduce the burden of repeated measurement, but the allocation of clusters and participants across treatment sequences requires careful consideration. We develop an optimal design framework for closed-cohort staircase trials with continuous outcomes, allowing cluster sizes to differ across sequences. Under a linear mixed model, we express treatment-effect precision through sequence-specific information contributions and establish convexity of the variance criterion in the cluster allocation proportions. An equivalence theorem characterizes optimal allocations, and reversal symmetry yields analytical solutions for several three- and four-sequence designs. We also consider joint optimization of sequence allocation and cluster sizes under a fixed participant total. Numerical investigations show that the value of unequal cluster sizes depends on the observation window and correlation structure: gains are often modest but can be substantial in selected settings. Comparisons with complete stepped-wedge designs demonstrate a trade-off between repeated measurements and additional clusters, with staircase designs achieving lower costs in some settings at comparable precision. Illustrations motivated by psychosocial cancer care and trial recruitment distinguish hypothetical design adaptations from the original studies. The results support choosing sequence allocation, cluster sizes, and observation windows jointly when repeated measurement is a major constraint.

Adjusting for baseline covariates has become standard practice in analyzing randomized clinical trials. In the low-dimensional setting, it is well understood that covariate adjustment through a parametric working model can sometimes be more efficient than the unadjusted difference-in-mean estimator. However, when the number of adjusted covariates is large relative to the sample size $n$, a naïve adjustment may introduce excessive bias, leading to invalid statistical inference. The current literature that tries to resolve this issue is either limited to linear working models or relies on sample splitting, which may raise concerns about the replicability of RCT analyses. In this paper, we devise a novel jackknife-based approach to covariate adjustment through generalized linear models (GLMs), which we term as JAckknife Score-based Adjustment (JASA), together with its calibrated version JASACal. By employing a nuanced jackknife strategy, JASA and JASACal avoid sample splitting and make full use of the data, while ensuring that the bias of JASA or JASACal is still negligible even when the number of adjusted covariates is large compared to $n$. JASA also encompasses state-of-the-art adjusted estimators through linear working models as a special case. Through extensive simulation experiments and a real data analysis, we demonstrate that JASA or JASACal can adjust for a much greater number of covariates than existing benchmarks. These empirical results also shed some new light on practical guidance for covariate adjustment with GLMs. Both JASA and JASACal have been incorporated into our R package HOIFCar available from CRAN. The package HOIFCar is developed to serve as a user-friendly option for covariate adjustment in RCTs, in particular when practitioners hope to adjust for a large number of covariates.

A switchback experiment alternates an entire system between treatment and control over time. It is especially useful when interactions between units can undermine standard unit-level experiments. Switchback experiments are commonly implemented on fixed temporal grids, which impose highly structured restrictions on when treatment can switch. We study a broader class of random-duration switchbacks, in which treatment and control alternate across runs whose durations are drawn from a common distribution. Under finite carryover, we show that the mean squared error of the Horvitz-Thompson estimator has a simple frequency-domain representation governed by how temporal outcome patterns align with design-induced imbalance and contamination patterns. This representation yields a tractable worst-case design criterion that can be computed directly from the duration distribution. Optimizing over even a simple two-rate family produces a run-age-dependent switching rule that reduces the asymptotic worst-case mean squared error by at least 32% relative to the standard independently randomized fixed-block switchback.

On evolving graphs, node classifiers must satisfy two key requirements: inductive generalization to newly arriving nodes under distribution shift and calibrated uncertainty for safety-sensitive applications. Standard graph neural networks (GNNs) are typically trained once and address neither requirement. We adapt the Bayesian last-layer (BLL) model by placing random last-layer parameters on top of a deterministic GNN encoder for uncertainty quantification. The categorical softmax likelihood required for classification breaks Gaussian conjugacy, so neither the training posterior nor the test-time streaming update has a closed-form solution. To address both challenges, we introduce a variational Bayesian last-layer (VBLL) objective that jointly trains the encoder and an approximate last-layer posterior by maximizing an evidence lower bound with a Monte Carlo expected log-likelihood. At test time, we freeze the encoder and apply an online Laplace update to the last-layer posterior. This update corresponds to a power-prior Bayesian model with exponential forgetting and a Kullback-Leibler anchor to the training posterior. Across five node-classification benchmarks under distribution shift, online GVBLL is the only method to achieve the best accuracy and negative log-likelihood on every dataset. It improves accuracy by up to 17 percentage points on Cora and 14 percentage points on ogbn-arxiv over the strongest non-GVBLL baseline, while remaining competitive in calibration with MC Dropout, Deep Ensembles, Temperature Scaling, and Gaussian-process classifiers.

Network analysis for multivariate time series is popular in many fields, from neuroscience to seismology. The inverse spectral density is a common choice for time series network analysis due to its representation of the frequency domain correlation between two variables after removing the best linear predictor of all other variables. In many applications, the goal is to study how these networks change across different conditions. For example, in neuroscience, one might be interested in how the brain connectivity network changes before and after stimulation. Towards this goal, we develop an inference framework based on a direct estimate of the difference in two high-dimensional inverse spectral densities. We develop a new Gaussian approximation error bound for any de-biased D-trace estimation procedure which is then leveraged to both inform optimal window sizes of Welch's estimators of the spectral density and establish asymptotic normality of our de-biased D-trace estimator. Moreover, we develop an efficient algorithm based on a generalized D-trace estimation procedure to overcome the computational complexity of high-dimensional inference. The method is illustrated on synthetic data experiments and on experiments with electroencephalography data.

Long records of rainfall contain critical information about the temporal organization of precipitation within a year, including wet and dry spells. Detecting when this temporal pattern changes is central to climate-impact assessment and environmental monitoring.A common data structure features observations from m locations over T years. For location i and year t, the record is a vector of daily measurements $X_{i,t} = (X_{i,t,1}, ..., X_{i,t,D})^\top$, with D = 365. The goal is to identify an unknown year $t_0$ when the temporal rainfall pattern changes across the region. Because changes in the intra-annual pattern affect the frequency-domain behavior of the daily sequence, we represent each yearly curve through a spectral feature capturing the strength of temporal oscillation at a fixed frequency.Classical changepoint methods often rely on parametric likelihoods or Gaussian approximations, which are hard to justify for skewed and heavy-tailed rainfall data. Furthermore, many procedures only return a point estimate. A confidence set for the change year is crucial in environmental applications to quantify temporal uncertainty and distinguish sharp transitions from weak evidence.This paper develops a distribution-free framework for detecting a synchronized structural change in yearly rainfall patterns across independently monitored locations. For each location-year pair, we compute a fixed-frequency nonparametric spectral density estimate. These features form a $T \times m$ data matrix. By exploiting the exchangeability of yearly spectral features before and after the true changepoint, alongside spatial independence, we construct finite-sample valid confidence sets for the unknown change year.

AI-based item generation and NLP-based prediction of item parameters are producing item banks that are substantially larger, sparser, and more frequently updated than conventional banks. Hierarchical Bayesian item response theory (IRT) is a natural calibration framework for such banks, but the common practice of refitting the entire accumulated response history at each update is costly and can exceed available memory. We describe consensus calibration, a divide-and-conquer procedure that calibrates each time period independently and reconstructs the pooled posterior in two layers. First, the posterior draws of each period are mapped to a common metric by a robust characteristic-curve linking (Haebara) that is solved separately for each draw, which propagates the uncertainty of the linking transformation into the linked posteriors. Second, the linked item posteriors are combined as a product of Gaussian densities from which the population prior contributed by each period is removed and a single prior—obtained by consensus across the per-period population posteriors—is reinstated. The correction targets the posterior dispersion, not only its location. As evidence for consensus calibration, we compare it to a pooled single-run analysis on a large operational assessment in terms of item-parameter recovery, an uncertainty-by-exposure diagnostic, and the ability distributions.

We develop a causal framework for matrix completion under missing not at random (MNAR) data. Drawing on synthetic controls from the econometric panel data literature, our approach relaxes two assumptions common in MNAR matrix completion: positivity and independence of observation indicators. Unlike traditional panel data models, which often require prescribed block-sparse geometries, our framework accommodates flexible, heterogeneous observation patterns through target-specific local information structures. We propose synthetic nearest neighbors (SNN), a local synthetic-controls-inspired estimator, and establish finite-sample entrywise error bounds and consistency for mean recovery under suitable conditions. We further derive asymptotic normality under heteroskedastic noise and develop feasible entrywise inference. To estimate entry-specific noise variances, we apply the same local principle to squared outcomes, obtaining consistency under bounded noise and asymptotic unbiasedness under general subgaussian noise. Simulation studies corroborate the theoretical findings across a range of missingness designs and observation patterns.