Papers

Gaussian mixtures of regressions are commonly implemented via a Gibbs sampler. This Markov chain Monte Carlo (MCMC) algorithm can be computationally burdensome because of the need to update discrete-valued latent component allocation parameters whose dimension increases as the sample size increases. In this article, we propose applying the method of composition to a Gaussian finite mixture model with a Normal-Inverse-Gamma (NIG) prior which allows one to write the posterior distribution as the product of conditional distributions. Namely, the conditional distribution of parameters given the data and mixture labels, times the marginal posterior of the mixture labels. The conditional distribution of parameters given the data and mixture labels, can be sampled from directly, instead of using MCMC. The expression of the marginal posterior of the mixture labels is known up to a proportionality constant and we adapt existing approaches in Bayesian selective inference to constrain the space of component labels to those arising from preliminary estimators, which alleviates a commonly encountered bottleneck. In simulation studies, we consider several settings and compare several versions of our constrained mixture of NIG models to two different MCMC-based strategies and demonstrate their use on natality data from the CDC.

We prove that the Kim–Milman flow map enjoys favorable stability properties with respect to variations in the target measure, provided that one of the target measures is sufficiently regular. Our results include stability in relative entropy, and more notably, Lipschitz stability in the $2$-Wasserstein distance up to a logarithmic factor. We complement our results with a general existence theorem for these maps for any target measure with finite second moment.

To address the computational challenges arising from large-scale longitudinal data, an optimal Poisson subsampling algorithm is proposed for quantile regression. The proposed method can substantially alleviate computational burden. Under some regularity conditions, we derive the asymptotic properties of the estimators from weighted quantile generalized estimating equations. For practical implementation, an efficient algorithm is proposed for parameter estimation. Furthermore, asymptotic theory is established for penalized weighted smooth quantile generalized estimating equations, and regularized parameter estimation is performed within the optimal Poisson subsampling framework. Both numerical simulations and a real data application demonstrate that the proposed optimal Poisson subsampling algorithm outperforms the uniform Poisson subsampling algorithm, and the regularized estimation exhibits satisfactory performance as well.

Neuron degeneration is the underlying mechanism for the development of many diseases. Quantifying the association between increasing levels of toxic exposure and progressive neuronal damage is a critical component of understanding this development. We investigate this association by analyzing a novel dataset of ordinal neuronal damage scores derived from a series of toxicological assays of C. elegans, including variables such as toxicant concentration, maternal treatment, and direct chemical exposure. We propose a computationally efficient parameter-constrained Bayesian ordinal regression that captures the monotonic association between neuron damage scores and corresponding treatments. Power analysis via simulation studies reinforces the advantages of our model over standard alternatives used in existing work by practitioners. Analysis of the novel C. elegans assays indicates that maternal toxicity increases susceptibility in progeny, with the offspring generation exhibiting amplified neuronal damage upon later-life rotenone exposure even under mild parental developmental treatment.

We present a framework for online and adaptive forecasting and hierarchical reconciliation using linear regression models. We begin by formalizing hierarchies using graphs, and motivated by their structure, formulate a multivariate linear model using the matrix normal distribution to characterize residuals. Parameter estimation is posed as a ridge regression problem and applied to hierarchical forecast reconciliation. The connections between ridge regression, Bayesian estimation and shrinkage for hierarchical reconciliation are discussed, and results for uncertainty quantification in parameters and forecasts are provided. Based on the ridge regression formulation, a recursive inference scheme inspired by recursive least squares is described. The algorithm is implemented in the PyOnlineForecast package. Finally, the proposed methodology is demonstrated on a case study for district heating load forecasting using a temporal hierarchy. Our results provide a reference for implementation of forecast reconciliation via multivariate linear models in an online setting. The case study furthermore highlights practical considerations of using temporal hierarchies in an online setting and demonstrates the usefulness of the proposed framework and implementation, both for district heating load forecasting and more generally for online hierarchical forecasting.

In this paper we prove some concentration inequalities for two types of error probabilities in the Empirical Risk Principle (ERP) in statistical learning, which provide a lower bound and an upper bound for the minimal risk (in terms of the minimal empirical risk) with non-asymptotic high confidence. The usual boundedness condition of the empirical risk function is relaxed to the Gaussian or exponential integrability condition. The confidence of the lower bound of the minimal risk is shown to be independent of the number of training parameters and the dimension of the input vectors, allowing one to detect the deficiency of a learning machine efficiently; and the confidence of the upper bound of the minimal risk is proved to be high provided that the sample size $n$ is much greater than the box dimension of the parameter set $\Theta$ in the Orlicz metric $d_{\psi_1}$ associated with the risk functions. Our work is based on Talagrand's concentration inequalities (the sharp versions by Bousquet and Klein-Rio), transport-entropy inequalities and the recent progress in the theory of empirical processes and statistical learning.

The objective of modern early oncology dose-finding is to identify an optimal biological dose (OBD), rather than simply the maximum tolerated dose. In basket trials, the dose-toxicity and dose-efficacy relationships may differ across biomarker or disease-defined subtrials, so a single common dose from pooled analysis may be suboptimal. We propose a flexible exchangeability-non-exchangeability (EXNEX) dose finding design (DF-EXNEX design) for subtrial-specific OBD selection in basket phase I/II trials with binary toxicity and continuous efficacy endpoints. Patient toxicity is modelled by a monotone logistic regression and efficacy by a quadratic dose-response curve. Robust borrowing is introduced through extended EXNEX mixture priors on the subtrial-specific curve parameters, allowing the strength of borrowing to adapt to the similarity of subtrials. Dose recommendation is based on an admissible set defined by posterior safety and futility rules, and an OBD-oriented utility function combining toxicity and efficacy on comparable scales. The operating characteristics were evaluated in a large-scale simulation study for the basket trial with four subtrials and five dose levels, and 70 scenarios covering all non-redundant combinations of true subtrial-specific OBD locations. Results showed that, compared with a no-borrowing NEX design, the DF-EXNEX design can increase the correct OBD selection for most scenarios while reducing overly toxic recommendation as final OBD. The improvement increased with subtrial similarity due to robust information borrowing, but a small number of mixed low/high OBD scenarios showed negative or near-zero gains, consistent with occasional over-borrowing towards intermediate doses. These results support robust borrowing for subtrial-specific OBD finding while highlighting the need to monitor borrowing behaviour when true OBDs are widely separated.

We study Transformer-type residual layers under cross-entropy training through a continuous-depth mean field control viewpoint. Depth is treated as time, layer parameters as controls, and the residual Transformer recursion as an explicit Euler scheme for a controlled hidden-state flow. For fixed controls, we prove an $O(\varepsilon)$ pathwise approximation of finite-depth trajectories by the continuous flow and combine this with high-probability sampling bounds for the empirical cross-entropy risk. We formulate the limiting population problem as a first-order transport control problem for the law of hidden states and derive a Pontryagin condition whose terminal adjoint contains the softmax residual. We also give finite-class and metric-entropy uniform estimates, compare optimal values, and discuss existence, stability, continuous-to-discrete recovery, initialization, and range estimates for continuous minimizers.

We study mirror flows generated by a convex quadratic loss and a general convex lower semicontinuous mirror potential. We show that, when initialized near the boundary of the domain of the mirror potential, their rescaled trajectories converge to a limiting mirror flow whose potential is the indicator function of the domain. In this limit, the primal variable minimizes the loss over a time-dependent hypothesis set: the subdifferential of the support function of the domain, evaluated at the dual variable. This characterization provides a general mechanism for incremental learning in mirror flows.

In high-dimensional data settings, dimensionality reduction or variable selection are key steps when using statistical learning techniques. Principal Covariate Regression-type methods aim to perform both dimensionality reduction and (regularized) regression steps in one analysis. However, existing PCovR methods cannot simultaneously select dimensionalities and estimate regularized coefficients, forcing researchers to make ad-hoc choices in the order of these steps. In this study, we propose a novel method called Principal Covariate Regression with Nuclear Norm Penalty (PcovRnnp) that allows simultaneous dimension reduction and estimation of regularized coefficients.

Health-related quality-of-life (HRQoL) outcomes are increasingly incorporated into oncology research to complement traditional survival endpoints by capturing patients' well-being over time. These outcomes are typically collected through multidimensional questionnaires yielding longitudinal ordinal data, and are often subject to dropout due to disease progression or death. In this context, joint models provide a well-established framework to account for the dependence between longitudinal HRQoL trajectories and time-to-event outcomes, but fully joint estimation rapidly becomes computationally prohibitive when multiple latent dimensions and random effects are involved. We propose a novel slope-corrected two-stage (SC2S) approach for the joint analysis of multivariate ordinal HRQoL data and survival outcomes within a multidimensional latent trait framework. The proposed approach propagates longitudinal information to the survival model through informative priors on the random effects, while additionally re-estimating longitudinal slope parameters. This strategy substantially reduces bias in both longitudinal and survival submodels while preserving much of the computational efficiency of two-stage procedures. Through simulation studies and an application to HRQoL data from patients with progressive glioblastoma, we show that the proposed method closely approximates fully joint Bayesian estimation while requiring notably less computation time.

Generalized linear models are central to actuarial modelling of binary risk, claim frequency, utilization, and cost-related outcomes. Yet fairness diagnostics often rely on linear-model intuitions, although GLM predictions are obtained by transporting a latent score through a nonlinear inverse link. We develop a moment-based decomposition framework for diagnosing group disparities in fitted GLM predictions. In an exact linear-Gaussian benchmark, the Wasserstein barycentric criterion for distributional demographic-parity violation reduces to a two-moment criterion and decomposes into direct mean, indirect mean, interaction, and structural components. For GLMs, we distinguish the empirical output-scale criterion $U_2(f)$, a within-group proxy $\widetilde U_2(f)$, and a leading decomposition $D_1(f)$. This leading term preserves the four linear channels and adds two curvature components induced by the inverse link: curvature coupling and curvature amplification. We derive explicit formulas for logistic, Poisson, and Tweedie specifications and illustrate the diagnostic on medical-expenditure survey data. The framework is not a legal test of discrimination, nor a full characterization of distributional parity outside the linear-Gaussian case. It is a tractable actuarial diagnostic for identifying whether fitted prediction disparities arise from explicit sensitive effects, proxy-mediated covariate profiles, covariance-structure differences, or nonlinear link effects.

Drawing on the observed regularity of typical trips, this study proposes a method for identifying commuting patterns without using the individual identifiers of public bike-share users. The approach combines an analysis of trip duration distributions, mixture models and logistic regression based on the posterior probabilities of belonging to the ___fast' and ___typical' journey categories.

The deployment of data-driven models in 6G wireless networks is increasingly challenged by frequent distribution shifts that degrade performance over time. Unsupervised Domain Adaptation (UDA) offers an alternative approach by adapting the trained model to a shifted domain without requiring labels. However, UDA pipelines are often more complex than single-task training due to additional modules and optimization procedures, raising a practical question: do the benefits of adaptation come at a higher energy cost, and how does this trade-off compare to retraining when labeling effort is also considered? In this work, we investigate the energy consumption of UDA and compare it to single task. We further propose a way to determine the minimum number of target domains for which UDA becomes more energy-efficient than retraining, taking into account the labeling cost. Our results aim to clarify when UDA should be preferred over classical train-from-scratch approaches from an energy and labeling-aware perspective.

Most tools for measuring political positions, manifesto coding, expert surveys, text-scaling models, were built and validated on Western party systems, and outside that setting they work poorly, and often not at all. This paper is an attempt at a method for those settings. It treats a large language model not as a measurement device but as a single, fallible rater in a panel, roughly the way an expert survey treats one expert: the value comes from pooling many judges rather than trusting any one of them. I describe the panel, an applicability rule that keeps a score of zero distinct from a blank, and a lens system that separates what an actor says from what it does. I report three results. First, holding a definition-free round fixed, adding written axis definitions moves scores by a mean of 1.8 points on a 21-point scale and tightens agreement between raters (mean absolute gap 2.81 to 2.50; r 0.81 to 0.89); they make two independent raters agree more closely, which an arbitrary steer would not. Second, across nine models from eight laboratories in two countries, Krippendorff's alpha is 0.86 on both an interval and an ordinal metric, and it stayed put as the panel grew from five raters to nine. That is reliability, the reproducibility of a reading, and not validity, its correctness. Third, where the panel does disagree, the disagreement is informative: the sharpest split, a full-scale divergence on an actor's stance toward its state's foundational order, points to a referent problem, and a blind triple-coding puts about two-thirds of it down to interpretation rather than error. I try to be plain about what the method can't do, including the human validation it still lacks, and I release the instrument and data in full. The worked example is the Middle East and North Africa, but I'd expect the method to carry to any region these standard tools leave out.

Functional data are often modeled through one likelihood-linked curve, while the scientific target is a larger state containing rates, accumulated quantities, boundary values, or nonlinear functionals of several linked levels. These targets require more than smoothing the observed curve: derivative uncertainty, cross-level covariance, and integration constants must be handled jointly. We introduce anchored Gaussian process differential ensembles, embedding an anchor \(f_0\) in a joint Gaussian state with its mean-square derivatives and repeated integrals. Integral levels add explicit Gaussian integration constants. This separates the anchor-induced covariance from finite-dimensional boundary uncertainty and clarifies why anchor-only observations do not identify independent integration constants. For stationary one-dimensional kernels, we compute the ensemble with a transformed Hilbert space Gaussian process approximation that applies derivative and integral operators to Laplacian–Dirichlet basis functions while retaining the integration-constant covariance exactly. We establish operator-level approximation bounds and conditional finite-grid posterior convergence. We introduce TARTARE, a target-aware calibration procedure for finite-rank differential ensemble approximations, to address derivative under-resolution by anchor-calibrated bases. In second-order simulations, derivative-aware calibration improves derivative posterior recovery relative to anchor-only calibration while preserving anchor and integral summaries. A motorcycle crash analysis illustrates coherent posterior inference on a coupled kinematic state and short-horizon turning-point functionals.

Decision-making in hierarchical systems requires probabilistic forecasts at all cross-sectional levels. Current hierarchical forecasting methods typically generate independent forecasts at each level and reconcile them post hoc to ensure coherence between upper and lower levels. Such post hoc corrections do not incorporate hierarchical structure or decision goals into the underlying parameter estimation. We propose a fully Bayesian hierarchical forecasting framework that shares information more effectively between and across levels than reconciliation alone. Our approach has the flexibility to softly penalise incoherence, subject to model specification, and to focus the global model and coherence update on hierarchical levels most relevant to decision outcomes. This yields parameter estimates that are focused towards the forecasting goals and capture the requirement for coherency, removing the need to estimate covariance matrices for multi-step forecasting horizons. We demonstrate improvements in predictive accuracy metrics on both simulated data and Australian domestic tourism forecasting.

We develop a comprehensive theory for regularized M-estimation in reproducing kernel Hilbert spaces. Under mild conditions on the loss we establish existence and measurability of the estimator, covering a wide range of convex and non-convex losses, including bounded robust losses. We further prove sharp rates of convergence with an explicit bias-variance decomposition governed by a novel complexity measure. We show that the variance is independent of misspecification, while the bias depends on a source condition parameter known in the learning literature. For tensor product Sobolev spaces we obtain new rates that connect to spaces of functions with dominating mixed smoothness, substantially extending existing results and explaining why these estimators circumvent the curse of dimensionality. Our methodology, combining elements from both functional analysis and empirical process theory, allows for an asymptotic linearisation of the objective function that avoids both closed-form solutions and global Lipschitz assumptions, and may be of independent interest. The estimators are implemented in C++ and theory is supported by numerical experiments.

Likelihood-free inference methods can perform Bayesian inference when evaluating the likelihood is impractical but simulating synthetic data from the model is feasible. Approximate Bayesian computation (ABC) is a well-established likelihood-free approach that constructs particle posterior approximations by evaluating the similarity between simulated and observed data using a distance function, which is used in rejection or weighting steps. Here we extend previous work on adaptive distance learning for ABC to misspecified time series, while also exploring applications in neural posterior estimation using prior-data fitted networks (NPE-PFN) with localization. The adaptation of the distance that we consider optimizes out-of-sample predictive performance using a scoring rule. We also establish a connection between linear pooling for forecast combination and our posterior estimation methods with randomized distances, showing that empirical estimation of pooling weights can be interpreted as another form of adaptive distance learning. For both ABC algorithms and NPE-PFN methods with localization, adaptive distance learning improves forecasting performance in simulated and real examples.

This paper studies kurtosis in multivariate normal variance-mean mixtures through its fourth-cumulant representation. We obtain an explicit expression for the fourth cumulant whose structure separates naturally into a rank-one directional component, a mixed direction–covariance component, and a covariance-pairing component induced by the mixing variable. This formulation shows that kurtosis in this class is not merely a directional tail phenomenon, but also reflects the interaction between mean variation, covariance structure, and stochastic mixing. We further derive the standardized fourth cumulant, relate it to Mardia's multivariate excess kurtosis, and study directional excess kurtosis through projection pursuit. Statistical applications are developed for cumulant-based diagnostics of multivariate non-Gaussianity, dominant-tail-direction analysis, and influential-tail-event detection. The practical relevance of the theoretical results is illustrated with simulated data and daily stock returns.

Nested sampling is a Monte Carlo algorithm for posterior estimation and Bayesian model comparison. It maintains a population of $K$ live points sampled from the prior, and at each iteration discards the lowest-likelihood point and replaces it with a new sample drawn from the prior restricted to exceed the discarded likelihood. Achieving this likelihood-restricted prior sampling efficiently and reliably is the central computational challenge. For low-to-moderate dimensional problems, MLFriends is a general and robust region-based approach that constructs a proposal region by bootstrap aggregation over the current live points and rejects proposals outside this region. We present a self-contained mathematical formulation of MLFriends and derive, under a homogeneous Binomial point process model for the live points, heuristic bounds on the expected fraction of the likelihood-restricted prior not covered by the proposal region. These bounds decay as $(\frac{1}{3}Km)^{-3/2}$, where $m$ is the number of bootstrap rounds, and are negligibly small for practical parameter choices. We show heuristically that the resulting bias in the marginal likelihood estimate is negligible compared to the inherent statistical variance of a nested sampling run. While a fully rigorous treatment remains an open problem, these results provide the first analytical characterisation of a fully specified and practically implementable nested sampling algorithm, without assuming an idealised or asymptotic sampling procedure.

Nested Sampling is a Monte Carlo algorithm enabling posterior estimation and Bayesian model comparison, and is especially robust in multi-modal posteriors. This is because nested sampling maintains a population of live points sampled from the entire prior. In each iteration, the population is advanced above a likelihood threshold, potentially discarding modes ruled out by the data. However, the Monte Carlo nature of point replenishment can also accidentally discard a mode. We draw a connection to the neutral Moran process in genetics, and quantify the occurrence probability of this failure mode of nested sampling with a simple symmetric random walk model on the live point occupancy. We find a simple rule for setting the minimum number of live points so that mode die-out is made unlikely.

Learning instability is a long-standing problem across machine learning, but it is especially acute in the overparameterized regime that defines modern deep learning: large models fine-tuned or trained on limited data traverse flat loss landscapes with many nearly-equivalent minima, and stochastic factors (initialization, data order, dropout, hardware non-determinism) can route optimization to very different solutions. The rise of large pretrained models (LPMs) makes the problem more urgent: training cost is high, downstream data is often small, and repeated runs for variance reduction are prohibitive. We introduce GRAIN (Group Aggregation via mIN-norm objective), a lightweight training algorithm that replaces the mean aggregation used in mini-batch optimization (both across mini-batches and within a mini-batch) with a min-norm convex combination of group-wise gradients. \mName guarantees a non-negative inner product between the aggregated update and every group gradient, resolving intra- and inner-batch gradient conflict, and retains an $\mathcal{O}(1/T)$ convergence rate comparable to SGD. Under mild smoothness and absolute-continuity assumptions, the min-norm solution differs almost surely from the arithmetic mean, which yields a uniform-stability bound for \mName strictly tighter than the standard bound for SGD. Empirically across generation, classification, and regression at LPM scale, \mName delivers consistent improvements in mean performance and reductions in run-to-run variance over a broad suite of tasks, with no extra training-time or storage cost beyond a single backward pass.

Density regression extends conventional parametric regression by allowing the entire distribution of the response to vary flexibly with covariates rather than just low-order moments. In the Bayesian setting, logistic Gaussian process (GP) priors have been widely used for density estimation and extend naturally to density regression. The prior can be centred on a base density model, with the nonparametric component providing an interpretable correction that is useful for model criticism. However, logistic GP density regression models have seen limited use, since they require computation of a normalizing constant for every observation, typically via numerical integration. We address this difficulty by proposing a generalized Bayesian approach using a loss function based on the Hyvärinen score. The Hyvärinen score depends only on derivatives of the log density with respect to the response, eliminating the need to compute normalizing constants. Since GP computations remain expensive, we also employ sparse inducing point approximations and variational inference to develop a scalable approach. We demonstrate the method on one simulated and two real datasets, including a German weather dataset with more than 150,000 observations.

Bayesian modelling workflows often consider multiple candidate models of varying complexity. Model selection is commonly used to navigate potential trade-offs between model complexity and generalisability to new data. We study when model selection is unnecessary or can even be harmful for predictive performance in finite data regimes and find that the need for selecting simpler models can depend on prior choice. We formalise predictively consistent priors, which keep prior predictive implications stable as model complexity increases. Across examples and numerical experiments, including adding covariates in linear and logistic regression, forward variable selection, and nonlinear modelling, flexible models with predictively consistent priors typically match or outperform selected simpler models in out-of-sample predictive performance. When selection helps, it can indicate poor joint prior implications, such as excessive prior mass on implausible predictive values. Based on our findings, we propose replacing the notion of sparsity or parsimony at the level of model components with specifying priors that remain sensible in predictive space as models become more complex.

Semiparametric efficiency theory provides the mathematical foundation for influence-function-based estimation, including one-step estimators, targeted minimum loss estimators, and many modern inferential methods used in causal inference and missing data analysis. Despite its widespread use, the theory is often presented through a collection of technical constructions whose geometric meaning remains opaque. As a result, influence functions are often derived and applied without an intuitive understanding of the principles connecting scores, tangent spaces, nuisance tangent spaces, and efficient influence functions. This tutorial develops a geometric exposition of semiparametric efficiency theory as a form of differential calculus on a space of probability distributions. Drawing systematic parallels with ordinary multivariable calculus, we show that paths of distributions play the role of curves, scores play the role of velocity vectors, influence functions play the role of gradients, and efficient influence functions arise as projected gradients. This perspective provides a unified explanation for several foundational questions, including why perturbation directions are represented by functions, why tangent spaces depend only on the statistical model whereas nuisance tangent spaces depend on the parameter of interest, and why efficient influence functions arise through orthogonal projection. The resulting framework offers a geometric perspective on semiparametric efficiency theory and influence-function-based inference.

Distribution shift between training and deployment is a pervasive challenge for modern AI systems. In many cases, the target marginals of covariates and response are known or specified through population-level observations, boundary conditions, properties of simulator configurations, or alignment-time distributional constraints. Such knowledge may provide valuable side information for regression estimation. We study this problem in the multivariate linear regression setting with a stable conditional mean $E[Y\mid X]$ across source and target, and identify the hybrid-loss estimator, which jointly incorporates both target marginals, as a benchmark target-aware estimator. Its direct computation, however, requires solving a coupled nonlinear optimization that is expensive at scale. Our main contribution is to develop and evaluate two computationally tractable alternatives: a constrained moment-matching estimator and a two-stage estimator that augments ordinary least squares with a calibration step. For all three estimators, we derive and compare closed-form asymptotic mean squared errors, yielding conditions under which the tractable alternatives match or closely approximate the hybrid benchmark, and regimes in which they do not. Monte Carlo experiments across three controlled shift regimes validate the theoretical results, investigate the accuracy-runtime tradeoffs among the three estimators, and translate into guidance on estimator choice. In particular, the two-stage estimator nearly matches the hybrid benchmark in the high signal-to-noise regime at essentially no additional cost, providing theoretical grounding for empirical observations in nonlinear settings.

Statistical matching combines partially overlapping datasets that share covariates $X$ but observe the target $Y$ and auxiliary variables $Z$ separately. Classical approaches typically invoke the conditional independence assumption (CIA), which makes the problem identifiable but fundamentally implies that the imported auxiliary variable provides no additional predictive power for $Y$ once $X$ is known. To capture this latent $Y$–$Z$ dependence, we propose a novel dependency-aware Schrödinger bridge for predictive statistical matching. Our approach couples the two separated databases by tilting the conservative CIA baseline with a transportation-based compatibility cost, recovering an informative joint distribution. The resulting statistical learning framework yields full probabilistic posterior rules for bidirectional imputation. Theoretically, we establish a sufficient condition under which the learned bridge strictly improves over the CIA baseline, alongside an exact joint recovery guarantee in the Gaussian setting under an appropriate cost. Across synthetic benchmarks and real-world datasets (CelebA and Adult), we demonstrate that our dependency-aware completion consistently improves downstream predictive utility, proving especially beneficial in settings like data recoding where the underlying population exhibits strong $Y$–$Z$ dependence.

The ability to communicate statistical results to domain experts and stakeholders is an important goal of the undergraduate statistics and data science curriculum. However, as large language models (LLMs) have become more accessible, a major concern is that students are offloading important cognitive tasks to generative AI. Using a corpus of over 1,600 undergraduate students' data analysis reports from 2021 to 2025, we show how students' writing style and verb usage have become more similar to that of LLMs. This shift is most pronounced in the first and fifth quintiles of students' reports, which roughly map onto the introduction and conclusion sections, respectively. At the same time, we demonstrate that students' writing style has become more similar to that of statistics experts with the addition of LLMs. We end by discussing the implications of our findings for statistics and data science educators. In particular, we propose alternative modes of assessment that still emphasize statistical thinking, such as targeted writing assignments for structuring a report introduction.

Likelihood is standard for Hawkes-process inference, while less computationally demanding methods have largely developed separately. We show that least squares, Takács–Fiksel, and related moment-based estimators form a single class of compensator-based estimating equations, with the likelihood score as the efficient benchmark. For fixed-dimensional multivariate Hawkes processes with compact memory, nonlinear positive links, and signed kernels allowing inhibition, every suitably regular predictable functional of a fixed lag window yields an unbiased estimating equation when integrated against $\text{d} N-\lambda\,\text{d} t$. Under common regularity, identification, and rank conditions, estimators based on every admissible finite library achieve uniform high-probability and pointwise almost-sure $\mathcal O(\sqrt{\log(T)/T})$ rates, asymptotic normality with Godambe covariance, and admit feasible two-step optimal weighting. A projection identity quantifies each library's exact efficiency loss as the score information outside its predictable span; a two-point bound shows the root-$T$ scale cannot be improved uniformly. Although compact memory localizes the intensity rather than the stationary law, exponential forgetting yields Bernstein-type concentration and transfers the theory to nonstationary starts after a logarithmic burn-in. Within this scope, the compensator class is exhaustive for finite-library comparisons: it contains the score, gives admissible libraries common guarantees, and quantifies their efficiency gaps exactly.