Papers

Mediation analysis plays an essential role in uncovering the mechanisms by which an exposure influences an outcome through intermediate pathways. While methodological advances for single-mediator settings are well established, rigorous tools for handling multiple, sequentially ordered mediators remain underdeveloped. Such settings are common in applications like longitudinal cohort studies, where exposures operate through complex chains of mediators over time. In this paper, we establish a general framework for sequentially ordered mediators that enables the identification and formal decomposition of the total effect into component path-specific effects. We also develop estimation procedures for mediation estimands with both continuous and categorical outcomes. Furthermore, we introduce a new testing strategy to conduct inference using a studentized statistic combined with data-splitting. This approach achieves valid Type I error control under the composite null across diverse data-generating mechanisms. Through extensive simulations and applications to two large-scale empirical studies, we demonstrate that the proposed methodology provides reliable estimation, valid inference, and improved power for discovering novel mediation pathways.

Hedonic price models are widely used to assess how environmental amenities affect property values, yet methodological guidance for estimating direct price effects remains sparse. We conduct an empirical Monte Carlo simulation to evaluate the performance of traditional and causal machine learning approaches for estimating the direct unmediated price effect of spatially delineated amenities on treated properties (DUET), a conservative lower-bound approximation for welfare changes with direct applications to benefit-cost analysis. Where previous simulations rely on parametric assumptions, we retain the actual data-generating process underlying over 1 million property transactions from upstate New York (1990–2024). By randomly assigning "treatment locations" across iterations we establish a "ground truth" that allows us to precisely measure estimation error. Our results demonstrate that generalized difference-in-differences (DID) regression consistently outperforms baseline DID and two-way fixed effects models across all scenarios. Causal Machine Learning (CML) methods, particularly causal forest DID, achieve comparable performance to generalized DID in most scenarios. In larger samples (above 3,000 treated) increasingly common in contemporary hedonic studies, CML approaches offer substantial advantages when properly specified. Based on empirical simulation results, we provide a set of method-specific best practice recommendations for both traditional regression and causal machine learning approaches.

Computer experiments with both quantitative and qualitative inputs have become common across various areas. However, constructing accurate and computationally efficient emulators for such experiments at large scales remains a significant challenge. We propose a novel, scalable framework for emulating computer experiments with mixed inputs. Our approach is based on a new covariance function integrating additive Gaussian Processes (GPs) to handle the mixed inputs, with Vecchia approximation for scalability. We demonstrate that methods for large-scale computer experiments can be effectively extended when paired with our proposed modeling framework.

The ASCCR (Attitude-Structure-Content-Communication-Relationship) framework was recently developed to teach collaboration skills to statisticians and data scientists. However, its effectiveness in real-world settings has not yet been systematically evaluated. To assess this, we evaluated novice undergraduate and graduate students' performances in initial collaboration meetings with real domain experts and compared them to an expert collaborator. Using video recordings, rubric scores, and domain expert feedback surveys, we found that novices performed surprisingly well compared to the expert. Specifically, novices scored nearly as well as the expert on the Attitude, Structure, and Relationship components of the ASCCR framework. Although novices did not initially perform as well on the Content or Communication aspects, they were able to close the gap. By the end of the collaboration projects, the novices had higher overall domain expert feedback scores than the expert. The primary implication of our study is that novices can become effective collaborators in a very short time. We discuss our findings' practical implications and provide recommendations for integrating the ASCCR framework into statistics and data science collaboration, consulting, and capstone courses.

Gradient boosted decision trees require a stopping rule to avoid overfitting. The standard rule monitors a validation loss and stops if the loss fails to improve for a fixed patience period. However, the patience parameter has no interpretable scale and validation losses can be noisy or implicitly defined by a user-specified gradient. We propose ScoreStop, a gradient-based early-stopping rule that casts the stopping decision at each iteration as a test of the null hypothesis that the current predictor is the population risk minimizer. We use a functional score test, computed on validation data, with a statistic that is scale-invariant in the update direction, with a known asymptotic distribution under the null. Because our test uses gradients rather than loss values, the same construction applies to implicit losses such as LambdaRank, and data-dependent losses such as Cox regression via influence functions. In synthetic experiments and real-data benchmarks, we show that ScoreStop is competitive with loss-based methods.

Visual and quantitative goodness-of-fit diagnostics are an important tool in the practitioner's toolbox. The need for convincing and reliable diagnostics is particularly clear when fitting extreme value regression models, which are used for extrapolation far beyond the observable range of the response variable, and often evaluated at unobserved covariate values. Despite this, few diagnostics have been developed for extreme value regression models, and those available often suffer in terms of interpretability or scalability on low-dimensional or non-Euclidean covariate domains, often encountered in modern applications. Moreover, existing methods tend to offer a global perspective on model fit; that is, they quantify goodness-of-fit across the entire dataset, without offering insight into regions of the covariate space where the model fit may be poor. We propose two novel visual diagnostics for extreme value regression models: the standardised tail plot and the normalised residual plot. By considering the asymptotic distribution of normalised exceedance probabilities, we show that uncertainty bounds for our plots are approximately independent of the sample size used in their construction. This allows us to propose visual diagnostics which can efficiently and consistently compare goodness-of-fit at both a global and regional level, despite varying sample sizes over regions of the covariate domain. Following a discussion of summary statistics for global and regional goodness-of-fit, we provide two applications of extreme value regression models that illustrate how our diagnostics can be used to perform model comparison (across thousands of candidate models) and provide actionable findings that support model design.

Latent state-space models are widely used to study partially observed dynamical systems, yet most formulations assume that process variability is independent of latent-state position. In many biological, behavioral, and physiological systems, however, variability may depend systematically on the underlying dynamical state, producing structured stochasticity that is not captured by constant-variance models. We introduce a state-coupled stochastic volatility framework in which latent process variance depends on displacement from a latent equilibrium. To estimate this relationship under partial observation, we develop a particle expectation-maximization procedure combining bootstrap particle filtering and backward trajectory smoothing. The model includes a coupling parameter, $\gamma$, that quantifies the strength of association between latent-state position and process variability. A large-scale simulation benchmark evaluated recovery and detection performance across varying coupling strengths, observation noise levels, trajectory lengths, and persistence regimes. The proposed framework consistently reduced recovery bias relative to an observed-state heteroskedastic proxy, with the largest improvements occurring under strong coupling. Recovery performance improved with increasing latent persistence, while detection performance remained competitive across a broad range of conditions and became increasingly advantageous as observation noise increased. Taken together, the results demonstrate that state-coupled volatility can be identified and estimated under partial observation when latent-state structure is explicitly modeled. The framework provides a practical methodological foundation for studying state-dependent variability and evaluating whether structured stochasticity contributes information about system dynamics beyond that contained in mean-state trajectories alone.

A fundamental goal in climate attribution is to estimate how forced climate change contributes to observed extreme weather events. The storyline attribution method compares an observed weather event, conditional on its atmospheric dynamic state (i.e., atmospheric circulation), in the current, 'factual' climate to an event with very similar circulation conditions in a hypothetical, 'counterfactual' climate. However, physical climate models cannot directly transfer these storyline counterfactuals across different climate forcing states. Statistical and machine learning techniques may overcome this limitation; yet, emulating circulation-conditional extreme events under different climate states is challenging. Here, we demonstrate distributional autoencoders (DAEs) as a versatile method for generating climate counterfactuals. They model the full distribution of spatially resolved European temperature fields conditional on the atmospheric circulation state and the mean global warming level. These distributions allow for deriving meaningful conditional probability ratios, which is a particular advantage of the DAE-based storyline approach. We train DAEs on fully coupled climate model simulations and we evaluate the modelled distributions across different factual and storyline-based counterfactual climate model simulations. In an illustrative case study, we revisit the 2003 European heatwave and we generate counterfactuals for a hypothetical `2003-like European heatwave' using ERA5 circulation, which we hypothesize to occur a quarter century (2028) and a half century (2053) after 2003. The conditional intensity would increase from 29.3 {\deg}C in 2003, to 30.3 {\deg}C and 32.1 {\deg}C in 2028 and 2053, respectively and conditional probability ratios would be 2.1 and 3.2 when compared to 2003.

Space-filling designs are commonly used in deterministic computer experiments. However, they are ineffective for factor screening, which makes them inefficient when only a small subset of input factors is influential to the output. Recently developed screening designs, such as MOFAT designs, are effective at identifying important factors but lack space-filling properties, limiting their usefulness for surrogate modeling. In this article, we propose a new class of screening designs that improves the space-fillingness while retaining their screening capability. Through several numerical examples, we demonstrate that the proposed designs offer clear advantages over existing designs.

In the last few years, AI-based models have become the centre of attention in weather forecasting due to their increasing accuracy and efficiency. Pioneering among weather services, ECMWF has developed its Artificial Intelligence Forecasting System (AIFS) model, which was first to provide data-driven ensemble forecasts in June 2024. Since July 2025, the AIFS ensemble model has been operational and runs in parallel with ECMWF's physics-based Integrated Forecasting System (IFS), which is considered the gold standard in weather prediction. The new AIFS model can generate forecasts ten times faster than the classical numerical weather prediction model, while consuming approximately a thousand times less energy. We present the results of our systematic assessment of the performance of the IFS and AIFS models by comparing the accuracy of raw and post-processed medium-range 10-m wind-speed ensemble forecasts generated operationally by the two models for the period between July and November 2025 for more than 9000 synoptic observation stations across the globe. The post-processed case involves the parametric ensemble model output statistics (EMOS) as well as the non-parametric quantile regression (QR) approach to correct any systematic inaccuracies in the raw forecasts. The predictive performance of raw IFS ensemble forecasts proves to be substantially superior to the skill of the raw AIFS predictions for all investigated forecast horizons. As expected, post-processing significantly improves the skill of both IFS and AIFS predictions, and, across most verification metrics, EMOS is superior to QR, especially for short lead times. Compared to the raw ensemble, the differences in skill between the matching IFS and AIFS predictions are substantially decreased by post-processing and are mostly significant at short lead times, when the IFS forecasts outperform their AIFS counterparts.

We introduce and study a new model of correlated uniform attachment (UA) trees, where correlation is sprinkled throughout the time evolution of the process. In this model, two UA trees are grown in parallel, and at each time step a new node is added to each tree, with an edge between it and a uniformly chosen existing vertex in the respective tree. The two choices of attachment are correlated: with probability $\alpha$, the edges attach to nodes with the same time label in both trees, and with probability $1-\alpha$, the choices are made independently. We study fundamental detection and estimation questions for this model, given two unlabeled trees. In our main result, we construct a consistent estimator of the correlation parameter $\alpha$, as the size of the trees goes to infinity. The construction of our statistic relies on two key ideas. First, we use Jordan centrality to identify subsets of vertices of each tree whose intersection has a sufficient number of common early vertices. The second idea is that, across multiple time scales, it is possible to approximately determine the labels of vertices that have attached to these early vertices, using the sizes of fringe subtrees. Our analysis includes novel quantitative bounds on the fraction of early vertices that remain central, which are of independent interest in the network archaeology literature.

Molecular dynamics (MD) is a key tool for simulating the dynamical behavior of atomic systems. However, MD is inherently serial, which makes it difficult to increase single-system throughput with concurrent compute. To address this, we introduce Langevin Speculative Dynamics (LSD), a distributed and model-agnostic speculative sampler for accelerating MD without adding relative error. Inspired by speculative methods in language and diffusion modeling, LSD uses a draft model to propose fast simulation steps and verifies them in parallel with a slower target model, applying a transport map from the draft to the target distribution. We extend speculative sampling to second-order Langevin dynamics, derive the achievable speedup as a function of physical parameters, show that LSD generalizes across different systems and draft-target combinations with a 3-9x speedup, and confirm theoretically and empirically that LSD samples trajectories from its target model distribution.

Two-arm phase II clinical trials often benefit from an interim analysis that allows early stopping for futility, but Bayesian calibration of such designs is usually based on computationally intensive Monte Carlo simulation. In this work, a simulation-free methodology is developed to obtain Bayesian optimal two-stage designs in two-arm phase II trials with binary endpoints using Bayes factors as the primary measure of evidence. Building on recent matrix-search methods for fixed-sample two-arm Bayes factor designs and earlier correction formulas for one-arm two-stage designs, the proposed approach derives exact expressions for the operating characteristics of a two-stage two-arm design with a single futility interim. Bayesian power and type-I error are obtained by correcting the corresponding fixed-sample quantities for trajectories that would have been removed by early stopping, yielding a fully numerical calibration procedure that avoids Monte Carlo error entirely. The resulting method searches over admissible interim and final sample sizes to identify the optimal design that satisfies target constraints on Bayesian power, type-I error, and the probability of compelling evidence in favour of the null hypothesis, while minimizing the expected sample size under the null hypothesis. The methodology is illustrated in realistic phase II settings, including a detailed re-analysis of the riociguat trial in systemic sclerosis. Overall, the approach extends simulation-free Bayes factor design methodology to the practically important setting of two-arm two-stage phase II trials and provides a transparent basis for Bayesian design calibration and sensitivity analysis.

Modern command, control, communications, computers, cyber, intelligence, surveillance, and reconnaissance (C5ISR) environments place substantial attentional demands on mission commanders. Failures in attention allocation in these high-risk settings can have severe operational consequences. This study investigates the efficacy of gaze-driven, attention-guided adaptive decision support tools, including visual-only and multimodal designs, in a high-fidelity simulated military command center. To characterize gaze and attentional dynamics during interaction with these tools, recurrence quantification analysis was applied to eye-tracking data. Stepwise regression using the Bayesian information criterion was then used to identify recurrence-based gaze metrics associated with performance. Results showed that the multimodal adaptive decision support tool was associated with significantly higher performance than the visual-only attention-guided tool. Average diagonal line length showed a negative linear association with performance, whereas entropy showed a positive linear association. Recurrence rate, determinism, and entropy also showed nonlinear quadratic relationships with performance. In particular, recurrence rate and determinism followed an inverted-U pattern consistent with the Yerkes-Dodson law. These findings suggest that effective performance in dynamic C5ISR contexts depends on a balance between structured and flexible visual scanning, and that recurrence-based gaze metrics can help characterize attentional dynamics during interaction with adaptive decision support systems.

We study sequential decision-making in partially observable environments against strategic, adaptive opponents, modeled as partially observable Markov games (POMGs). The central challenge is to learn latent dynamics from partial observations while facing an adversary whose behavior depends on the learner's strategy, making standard regret notions inadequate. We prove that an epoch-based optimistic maximum-likelihood algorithm achieves ${\tilde{O}}(\sqrt{T})$ policy regret for fixed problem parameters, with explicit dependence on the horizon, adversary memory, confidence radius, and the aggregate Eluder dimension of the observable-operator error classes. The algorithm deploys one policy per epoch, with geometrically capped epoch lengths, confidence sets built cumulatively from past data, and a statistical termination test that ends an epoch as soon as the data refute the deployed optimistic model. We also prove a matching lower bound, and extend the framework to horizon-adaptive guarantees and geometrically fading adversary memory.

Bayesian optimization (BO) is a popular and effective approach for tuning expensive, noisy experiments, but requires the formulation of an explicit objective function. Preferential BO (PBO) removes this requirement by learning from pairwise human feedback, yet existing methods struggle to efficiently optimize beyond low- and medium-dimensional problems due to their global search approaches. We address this limitation by developing a family of local PBO methods that transfer key ideas from high-dimensional BO to the preferential setting. In particular, we introduce local PBO methods which adapt trust-region and derivative-informed local search to pairwise preference feedback, where the latter exploits first- and second-order derivatives of the Laplace-approximated GP posterior. Our benchmark on GP sample paths, standard optimization benchmark functions, and policy-search tasks shows that local PBO methods are especially effective in high-dimensional and complex landscapes with steep optima. Compared with global preference-based baselines, they can substantially reduce cumulative regret, making them particularly useful for real-world preference-based optimization tasks such as policy search.

Many machine learning problems, including similarity learning, ranking, and clustering, rely on empirical pairwise loss functions whose quadratic computational cost quickly becomes prohibitive at scale. We demonstrate how a frugal approach that retains only a fraction of the available information on pairs can achieve estimation or optimization performance comparable to that obtained by using all pairs, by leveraging survey sampling techniques. A central finding, supported by both theory and experiments, is that such sampling plans must target pairs directly rather than individual observations. In particular, for pairwise losses between high-dimensional vectors such as embeddings in vision or graph learning, assigning higher inclusion probabilities to informative pairs using suitable auxiliary information yields performance close to full pairwise evaluation, providing a principled and theoretically grounded trade-off between accuracy and computational cost.

Paired-comparison data are naturally represented by tournaments, where transitivity corresponds to the existence of a global ranking consistent with all pairwise outcomes. Accordingly, the classical Kendall-Smith coefficient of consistency measures deviations from transitivity in a tournament by counting the number of circular triads (directed $3$-cycles). In this paper, we characterize the fluctuations of the number of circular triads in inhomogeneous random tournaments and develop an inferential framework for the consistency coefficient. Specifically, we consider the $W$-random tournament model, where the comparison probabilities are determined by a tournamenton $W$, the analogue of a graphon in the tournament setting. We show that, for a $W$-random tournament on $n$ vertices, the number of circular triads exhibits three different fluctuation regimes, determined by suitable notions of regularity and uniformity of $W$. We further develop a novel tournamenton multiplier bootstrap that consistently approximates the limiting distribution of the circular-triad count in the relevant asymptotic regime. Combining this with procedures for testing regularity and uniformity, we design an algorithm for constructing confidence intervals for the consistency coefficient that is asymptotically valid for all tournamentons. We also obtain structural characterizations of tournamentons for which the limiting distribution of the number of circular triads exhibits specific degeneracies. These results can also be viewed through the lens of tournament quasirandomness and may be of independent interest.

When it comes to estimating an unknown spectral density as simply and reliably as possible, parametric spectral density estimation using AR models and order selection via AIC is the method of choice. In contrast, no standard method has yet emerged for automatic nonparametric spectral density estimation, and there seems to be little willingness to weigh the advantages and disadvantages of different risk functions and the various methods for estimating them on a case-by-case basis, particularly because it is unclear whether the effort is even worthwhile without concrete prior information about the unknown spectral density. As a result, subjective visual methods are still widely used in practice to determine the appropriate smoothing parameter for a nonparametric estimation. This article aims to encourage the increased use of objective automatic methods by presenting evidence that using what is arguably the simplest and most straightforward frequency-domain version of the AIC for the automatic determination of an appropriate bandwidth enables results that are comparable to those obtained using the standard parametric approach. This evidence is based on both real-world time series and synthetic time series with spectral densities of varying complexity.

The deterministic incompressible Navier-Stokes equations are physically incomplete: any viscous fluid at finite temperature must exhibit thermal fluctuations whose form is dictated by the fluctuation-dissipation relation. We formulate the stochastic Navier-Stokes equations with the kinematically selected deformation Laplacian on compact Riemannian manifolds with strictly negative Ricci curvature. The fluctuation-dissipation relation, derived from a topological (Poincaré lemma) argument, uniquely determines the noise from the viscous operator. For the spectrally truncated system, we prove that the unique stationary distribution is the Gibbs measure (Gaussian in the mode amplitudes, because the nonlinear convective terms preserve energy), and that convergence to equilibrium is exponentially fast with rate at least $2\nu\lambda_\Def$, where $\nu$ is the kinematic viscosity and $\lambda_\Def$ is the spectral gap of the deformation Laplacian. The spectral gap satisfies $\lambda_\Def \geq \kappa^2$ when $\Ric \leq -\kappa^2 g$, and is independent of the volume of the domain. On flat space, the analogous thermalisation rate vanishes in the infinite-volume limit. The equilibrium velocity-velocity correlation function decays exponentially in geodesic distance, in contrast to the algebraic decay on flat space. These results provide a rigorous statistical-mechanical foundation for viscous fluids on negatively curved manifolds and illustrate how the geometry of the domain controls not only the deterministic dynamics but also the approach to thermal equilibrium.

Shapley values are a principled attribution measure widely used in interpretable machine learning, but their exact computation scales exponentially with the number of players, motivating a wide range of approximation methods based on value function evaluations of sampled coalitions. This raises the question of whether approximation accuracy can be improved by adaptively selecting coalitions for evaluation based on previous evaluations. This is particularly relevant in settings where the value function is costly and the number of evaluations is severely limited, such as retraining-based feature importance, data valuation, and hyperparameter importance. For this purpose, we propose ShaplEIG, a Bayesian experimental design approach that approximates the expensive value function using a Gaussian process surrogate and adaptively selects coalitions based on their expected information gain about the Shapley values. By the linearity of the Shapley values in the value function, we show that the expected information gain is available in closed form. Furthermore, we propose an efficient computation scheme that reduces the complexity from exponential to polynomial in the number of players via elementary symmetric polynomials. In extensive experiments across diverse costly applications, our method consistently improves sample efficiency in the low-budget regime over state-of-the-art baselines.

This paper clarifies the identifying assumptions underlying causal inference based on treatment changes rather than levels, and their relationship to conventional identification strategies. We characterize two structural models, with non-nested assumptions, under which treatment-change identification is valid conditional on observed covariates by differencing out time-constant confounders that are additive in the treatment equation. The assumptions underlying treatment changes are generally not nested with those of methods relying on treatment levels, such as selection-on-observables strategies that control for past outcomes, treatments, and covariates, or difference-in-differences approaches that difference outcomes rather than treatments over time. We show, however, that under a random-walk restriction on the treatment process, exploiting treatment changes for identification is equivalent to using treatment levels given lagged treatment. This and other equivalence results motivate overidentification tests based on methods considering treatment levels and changes. Under an alternative model that does not assume a random walk but instead rules out dynamic treatment effects (among other conditions), treatment changes can still be used as an instrument to identify a treatment effect that is constant given covariates. However, without random walk, different identification strategies are generally not nested. In partially linear models, the non-nesting results carry a double robustness implication for two-way fixed effects regression that differences both the outcome and the treatment over time, which under certain conditions remains consistent if either the treatment-change assumption or the parallel-trends assumption holds. We characterize the causal models consistent with each method, run simulations for illustration, and present an empirical application to cigarette demand.

Temporal systems often exhibit non-stationary behaviour, such as seasonal climate variation or glucose fluctuations in patients with type-1 diabetes. One way to model non-stationarity is through discrete latent regimes, i.e., stationary segments of time. Such systems induce a Markov Switching Model (MSM), a class of Hidden Markov Models with autoregressive dependencies among latent regimes and observed variables. Identifying latent regimes is challenging in the presence of frequent regime switches and nonlinear and non-Gaussian dynamics, particularly when there are instantaneous effects between the variables, e.g., due to slow rates of measurements. In this work, we establish the identifiability of both latent regimes and regime-dependent causal structures under temporal regime dependencies, nonlinear lagged and instantaneous effects, and independent noise from the exponential family. Our identifiability theory subsumes non-temporal mixtures of causal models. Furthermore, we introduce FlowMSM, a regime detection framework that can be paired with any stationary causal discovery method to recover regime-dependent causal structures. Experiments on synthetic benchmarks and a financial economics dataset demonstrate the effectiveness of our approach to detect latent regimes and discover causal structures from non-stationary time series.

Predicting whether an individual with Alzheimer's disease will experience mild or severe disease progression is essential for personalized treatment. Typically, practitioners seek to predict the distribution of a discrete disease score, conditional on an individual's current MRI volume and their historical disease trajectory. Classical statistical regression models and single-task neural networks are not well-suited for this purpose because fitting separate models is infeasible (since each individual typically has few observations), while ignoring individual-level correlation leads to poor generalization. Meta-learning, in contrast, provides a natural avenue to dynamically predict distributions without retraining and model nonlinear relationships between the outcome and covariates. Motivated by this, we propose a Bayesian meta-learner that is trained on multiple individuals but tailors the predictive disease score distribution to each individual's historical data. Our model predicts on unseen individuals without retraining, scales linearly with the number of historical observations, and is guaranteed to be less overconfident when predicting long-term disease scores compared to its deterministic counterpart. On real-world data from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database, our model achieves performance competitive with both single-task models and deterministic meta-learners, while substantially improving performance when predicting long-term disease progression.

Estimating the generative mechanism of large-scale networks is a fundamental challenge in statistical machine learning. It requires the identification of the latent connectivity structure, which is in general an NP-hard combinatorial problem due to the absence of canonical node labels. We address this challenge by allowing for probabilistic couplings, thereby relaxing the assignment problem. Our estimation framework can be formulated as a semi-relaxed Gromov-Wasserstein objective and provides a low-dimensional representation of the generative structure. We solve this via a block-coordinate conditional gradient algorithm. Despite the relaxation, the resulting solution is typically deterministic: in fact, we show that the optimality gap between the relaxed solution and the deterministic assignment vanishes at rate $O(1/n)$, where $n$ is the number of nodes. This allows for tractable recovery of the underlying model and enables rigorous statistical analysis: we establish consistency and minimax-optimal convergence rates for both stochastic block models and Holder-smooth graphons. Our implementation scales efficiently with $n$, as demonstrated on both synthetic and real-world datasets.

Multinomial count data, such as microbial composition profiles derived from sequencing studies, frequently contain anomalous observations that distort parameter estimates. The Dirichlet-multinomial (DM) distribution is widely used in this setting but remains sensitive to such contamination. We propose the contaminated Dirichlet-multinomial (CDM) distribution, a two-component mixture in which the regular data come from a DM component with a lower dispersion and the irregular data come from a DM component with an inflated dispersion parameter. This construction accommodates anomalies without requiring their removal, and yields a natural rule for anomaly detection via posterior probabilities. Through sensitivity analyses involving both single-point anomalies and background noise, we demonstrate that the CDM distribution effectively downweights the influence of anomalous observations on the parameter estimates. The model is applied to gut microbiome data from a colorectal carcinogenesis study, where it consistently outperforms the DM distribution across all information criteria and identifies biologically plausible anomaly proportions in both the healthy and carcinoma subsets.

We investigate support thresholds for fully smeary and directionally smeary absolutely continuous probability measures on the sphere \(\mathbb{S}^m\). The motivation is inferential: smeariness is caused by degeneracy of the Hessian of the Fréchet function, and such degeneracy can invalidate the classical central limit theorem (CLT) for Fréchet means and the corresponding Wald-type \(\chi^2\) inference. For rotationally symmetric densities, we show that full and directional smeariness are equivalent. The Hessian and fourth-order terms are governed by two explicit geometry-dependent radii \(R_m<S_m\). In dimensions \(m=2,3\), rotationally symmetric smeariness cannot occur. For \(m\ge4\), support contained in the geodesic ball of radius \(S_m\), centered at the Fréchet mean, rules out smeariness; conversely, for every \(\varepsilon>0\) with \(S_m+\varepsilon<\pi\), we construct examples of rotationally symmetric \(2\)-smeary densities supported in the ball of radius \(S_m+\varepsilon\). For general densities, closed hemispherical support rules out both full and directional smeariness. Support contained in the closed ball of radius \(S_m\) rules out full smeariness, while we construct explicit, directionally \(2\)-smeary examples supported in balls of radius \(\pi/2+\varepsilon\). As a byproduct, the explicit Hessian formulas in this paper also provide a practical diagnostic for detecting proximity to the Hessian-degenerate, non-classical regime.

This article describes the design of a neutral comparison study in the context of empirical studies where the interest is in learning the functional relationship between a continuous errorprone exposure variable and a binary outcome. The performance of combinations of measurement error correction methods and flexible regression modeling techniques was compared using a simulation study. The project involved four independent teams, one devoted to data generation and evaluation, the other three to specific methods of measurement error correction (Simulation-Extrapolation, Regression-Calibration and Multiple imputation, Bayesian method). The study was conducted in three successive stages. In Stage 1, the first team simulated five datasets differing only by the true exposure-outcome functional form and distribution of true exposure. Furthermore, the implementation of flexible modeling methods (B-splines, P-splines, and fractional polynomials) was standardized. The three methods teams, blinded to the underlying data generation process, created the codes to implement their methods, and provided their results to the first team who evaluated them. These codes were then used by this team in the next Stages of the project. In Stage 2, the team simulated 150 additional datasets where other design parameters varied while using the same five exposureoutcome functions. Stage 3 consisted of simulating independent replications of each of the 150 scenarios considered in Stage 2 to quantify the sampling variance of the estimates. This work emphasizes the relevance of neutral comparison studies to fairly evaluate statistical methods aimed at addressing a complex analytical challenge, and demonstrates their feasibility through a large collaborative project.

In this work we provide analytical and closed-form expressions for the exact computation of the score and the observed Fisher information matrix in a Gaussian random walk observed through Gaussian noise. Our method is based on the Oakes' identity and, as for the computation of the log-likelihood, its complexity in time is linear in the length of the sequence with the forward-backward (or Baum-Welch) algorithm. We illustrate the method over various simulation studies and provide parameter estimates computed with the Newton-Raphson algorithm along with confidence intervals.

Probabilistic time series forecasting has attracted increasing attention in financial applications due to the need to quantify risk and uncertainty in future observations. We propose ProbRes, a post-hoc probabilistic calibration method that explicitly learns and incorporates volatility dynamics into probabilistic forecasting, enabling effective handling of heteroskedastic data. During training, ProbRes employs two architecture-agnostic modules to separately model the conditional mean and conditional volatility. At the inference stage, it generates predictive distributions by resampling normalized residuals. ProbRes is applicable to both univariate and multivariate time series and remains robust under a wide range of error distributions, including non-Gaussian innovations with conditional heteroskedasticity. Theoretical results demonstrate ProbRes's validity and experiments on both synthetic and real-world datasets show that ProbRes accurately captures predictive distributions and produces well-calibrated prediction intervals.