Papers

High-dimensional multi-way (tensor) data pose significant challenges for covariance estimation due to the curse of dimensionality. We introduce a unified framework for scalable estimation of tensor covariances based on a Kronecker-structured sparse inverse Cholesky (KSIC) projection. Our approach is grounded in the geometry of information projection, defining the estimator as the moment-matching projection of a target distribution onto a manifold characterized by sparse, Kronecker-factored inverse Cholesky factors. By leveraging physical or data-driven nearest-neighbor sparsity, KSIC provides a geometry-aware representation that is both statistically interpretable and computationally efficient. Our framework integrates two estimation regimes: a nonparametric estimator that projects the empirical covariance directly onto the manifold, utilizing the KSIC structure to implicitly regularize rank-deficient data; and a parametric estimator that fits generative covariance models (e.g., Matérn) by maximizing the likelihood of their KSIC projections, formulated as a nested double forward Kullback-Leibler minimization. Theoretically, we establish the conditions for the existence of the KSIC projection and finite-sample concentration rates for the nonparametric regime, proving that the KSIC estimator gainfully exploits cross-mode information and is robust to data scarcity. Numerical experiments demonstrate that the proposed KSIC estimators achieve state-of-the-art accuracy and scalability, particularly in settings with high dimensionality and limited sample sizes. We apply KSIC to spatiotemporal temperature anomalies and functional MRI data, demonstrating its broad applicability across diverse multi-way data domains.

Modeling multivariate distributions with nonlinear dependence, multimodality, and tractable analytical structure for downstream applications is a central challenge in uncertainty quantification. Sliced Normal (SN) distributions were introduced in prior works at the National Aeronautics and Space Administration (NASA) to address this need by representing densities through polynomial feature maps. This construction provides a compact algebraic alternative to more opaque generative models, while retaining the ability to capture nonlinear parameter dependencies and multi-modal behavior. In this paper, we build on the SN framework and develop several improvements that make the approach more reliable and scalable. First, we reformulate SN parameter estimation as a convex optimization problem over a positive semidefinite matrix, replacing the original nonconvex likelihood search with a formulation amenable to standard optimization tools. Second, we clarify the expressive power of the SN class by connecting polynomial log-density modeling to a Stone–Weierstrass-type universal approximation argument on compact domains. Third, we propose a high-dimensional fitting procedure that partitions variables into approximately independent groups, fits SN models within each subgroup, and then assembles the subgroup models through a cross-block completion step to recover residual dependence. We demonstrate the resulting SN modeling pipeline on NASA loss-of-control flight data, where the method captures nonlinear dependence patterns in both low-dimensional slices and a higher-dimensional block-assembled model.

We propose a novel regularity-informed filtering framework for data assimilation in the context of hyperbolic conservation laws and other time-dependent partial differential equations. We focus on systems whose states exhibit steep gradients and jump discontinuities. While filtering is widely used to improve numerical simulations by incorporating observational data, traditional filtering methods lack awareness of the spatial regularity of states produced in these systems. As a result, data assimilation often produces unphysical state estimates, introducing spurious oscillations in smooth regions and smearing sharp features. To address this limitation, we introduce a filtering framework that incorporates edge-preserving regularization into the filter's analysis step; this framework balances simulation forecasts, observational data, and structural prior knowledge. We formalize this approach using the ensemble Kalman filter (EnKF) and a class of hierarchical generalized sparse Bayesian learning (GSBL) priors, which adaptively infer spatially varying hyperparameters to promote non-oscillatory behavior in smooth regions while preserving discontinuities. We demonstrate the effectiveness of the resulting GSBL-EnKF method on challenging benchmark problems governed by hyperbolic conservation laws. Our results show that enforcing regularity in the analysis step yields sharper, less oscillatory state estimates and lower errors of the ensemble mean. This sometimes comes at the cost of ensemble spread, which we quantify and discuss.

Synthetic tabular data are often judged by realism, privacy, or downstream-task scores. Those scores do not answer whether a proposed release is supported for a named use, population, and threat model. We introduce SynthGuard-ReleaseBench, an audit framework that locks the use, candidate panel, tolerances, and audit schedule before evaluation. It compares real-trained and synthetic-trained workflows on protected data, gives simultaneous finite-sample bounds for bounded loss gaps, requires controls, and keeps utility, empirical privacy risk, mechanism claims, and human release authority separate. Across four American Community Survey studies, five non-ACS records, two chronological diagnostics, and a sealed prototype, the benchmark retains favorable, unfavorable, and excluded outcomes. Transparent baselines pass some locked audits; compact learned models fail under the declared budgets; a health-table case is excluded because its negative control passes. A post-audit scaling arm, repeated across three generation seeds, shows the same locked criterion admitting those learned models once they are fit on enough data while still rejecting a dependence-destroying control at every size, so the criterion discriminates rather than merely rejects; the same repetition withdraws a finer single-seed ordering. The theory adds a pre-audit sample-size rule, variance-adaptive and anytime-valid certificates that tighten the bound two to ten times on the same locked evidence, a temporal certificate for time-ordered audits, and two lower bounds: ordinary bounded queries reconstruct a protected audit once the query budget reaches its size, and the panel-size correction is necessary rather than conservative. The contribution is a reproducible workflow for use-specific release evidence, not a claim that any generator is private, safe, or deployment-ready.

Modern generative models increasingly produce distribution-valued outputs, such as predicted cellular responses to genetic perturbations in single-cell genomics. While these models provide valuable auxiliary information, they are inherently imperfect, creating a need for statistical methods that leverage their predictions without relying on their correctness. We propose generation-powered inference (GPI), a general framework for improving inference on distribution-valued parameters using auxiliary generative models. Focusing on Wasserstein barycenters and related distributional functionals, we introduce a function-valued bridge representation that transforms inference in the nonlinear Wasserstein space into estimation of a mean function in a Hilbert space, enabling an augmented estimation framework analogous to prediction-powered inference. We develop a family of GPI estimators with optimal information borrowing, establish consistency, asymptotic normality, and simultaneous confidence bands, and derive valid inference for linear functionals and Wasserstein distances. Simulation studies demonstrate efficiency gains over labeled-data-only methods and robust performance under generative model misspecification. We illustrate the proposed framework using a Perturb-seq study of K562 cells, where synthetic perturbation responses generated by the State foundation model are used to improve inference for pathway-level consensus gene expression distributions associated with perturbations of the 40S ribosome module.

Network datasets in modern applications often involve multiple types of interactions occurring over a shared set of individuals. Characterizing the generating mechanisms of these interactions can be enhanced by joint modelling, as shared vertices allow layers to help explain the structure of other layers. We model multiplex observations using graph limits, called a scaled set of graphons, and develop a nonparametric joint estimator based on blockmodel approximations, termed the multi-network histogram. This nonparametric framework captures each layer's varying sparsity and connection structure, accounting for heterogeneity via shared latent variables across all layers. We establish the theoretical properties of the multi-network histogram, providing an upper bound for the weighted mean integrated squared error and deriving the optimal bandwidth that minimizes this error. By leveraging information across layers, this joint modelling achieves a reduction in error and a smaller optimal bandwidth, which enables high-resolution estimation even in sparser layers. Its usefulness is demonstrated through simulation studies and an application to socioeconomic networks in an Indian village.

This study investigates the methodological and theoretical properties of session handover in applications that use large language models. A task may continue in a new session when the context reaches the model's input limit, when the application restarts, or when another agent is asked to finish the task. The application must then decide which information from the earlier session to pass on. We formulate handover as the transfer of a task-relative in-context learning (ICL) state and distinguish exact recovery of earlier material from preservation of the target distribution. Under an exogeneity condition, predictive equivalence characterizes the coarsest deterministic sufficient handover and gives a fixed-length bit requirement. The analysis isolates the effects of the memory constraint, the writer, and the continuation procedure, and quantifies the cost of writing before the realized downstream query is known. We propose a three-part record that stores decisions and constraints exactly, uses task-justified statistics for repeated evidence, and retains original observations whose effect is not preserved by those statistics. Gaussian linear regression gives an exact finite-dimensional handover and finite-bit perturbation bounds, while nonparametric regression gives upper and lower bounds that relate memory to squared prediction error. These results provide a theory and method for deciding what a handover must retain and how its memory requirement depends on the continuation task.

Whittle-Matérn (WM) Gaussian random fields (GRFs) are defined as solutions of stochastic partial differential equations (SPDEs) and provide a natural analog of Matérn GRFs on non-Euclidean geometry where the Matérn covariance function is not valid. In particular, WM GRFs on metric graphs have been an active area of research motivated by road and river networks where spatial dependence is more naturally described by intrinsic distances in the network than by Euclidean distances. This family of GRFs is controlled by three parameters relating to marginal variance, spatial range, and smoothness, but can be extended to so-called generalized WM GRFs through spatially varying coefficients in the SPDE. Recent work has considered the use of spatially varying covariates, but the full possibilities of flexibility have not been considered. In this work, we introduce latent GRFs that describe the spatially varying coefficients of the SPDE. This flexible model is compared to less flexible models in a simulation study evaluating both the ability to estimate the covariance structure and predictive ability. An important focus is the number of observations and replications necessary to reliably recover the covariance structure. We find that the flexible model improves over less flexible models in the presence of sufficient data. We also demonstrate practical applicability on traffic counts in a part of Madrid, and observe major differences between in-sample and out-of-sample predictive abilities of the models compared.

We present new inference tools for change point detection in high-dimensional time series. We discuss two distinct statistical applications: First, sequential change point testing in an incoming data-stream. Second, retrospective localization of multiple changes, with confidence intervals at a globally controlled error level. Test statistics are built on the maximum norm to generate power against sparse and asynchronous changes. Both problems are tackled by related multiscale statistics that search for changes in the data at many different levels of resolution. For fixed dimension, our statistical approaches can be validated using traditional Hölderian invariance principles. In this paper, we present the high-dimensional analogue: Hölder-Gauss-approximations, which can be (roughly) interpreted as the Gaussian approximation for a Hölder-norm of the high-dimensional partial sum process. Such approximations are of interest beyond change point detection and can be used for other problems such as for stationarity testing in high dimensions. We evaluate finite-sample performance in a simulation study and give an application to air contamination due to wildfires in California, which occurs asynchronously across a panel of measuring stations.

Ensuring food security across all regions of a country requires continuous monitoring, yet household surveys often leave significant spatio-temporal gaps due to resource constraints and operational priorities. In this paper, we propose a spatio-temporal additive Gaussian process model to estimate sub-national food security time series by regions. To address the computational cost of Gaussian process models, we exploit Kronecker structure of the spatio-temporal covariance matrix for scalable inference. We evaluate the proposed approach on food security survey data from Nigeria and Chad comparing it against other statistical and machine learning models and show how our proposal achieves better accuracy while retaining reliable uncertainty, especially when covariates are informative. We further apply the model to generate estimates for Nigerian states not covered by the survey, demonstrating its operational value for filling geographic gaps in food security monitoring.

Network data obtained through link-tracing designs, such as snowball sampling, are collected through a mechanism that depends on the very structure the analysis seeks to estimate. Ignoring this dependence and treating the observed sample as though it were itself a complete network can lead to substantially biased inference. While the resulting selection problem is intractable in general, we show that it admits an exact solution for $r$-wave snowball samples, with full-neighbourhood recruitment, drawn from an Erdős–Rényi population. We derive the exact likelihood of such a sample and show that it defines a curved exponential family in the edge probability $\pi$, with a low-dimensional sufficient statistic. Building on this result, we obtain the maximum likelihood estimator of $\pi$ that correctly accounts for the sampling design and, as a function of the minimal sufficient statistic, makes full use of the information in the sample. Simulation studies show that this correction substantially reduces bias relative to the naive estimator, remaining effectively unbiased even when the sample covers as little as 0.1% of the network. We further construct valid confidence intervals for $\pi$ by inverting a test built on the exact sampling distribution, approximated via Monte Carlo simulation. Simulation studies confirm that these confidence intervals attain the nominal coverage level within Monte Carlo error across a range of edge probabilities and numbers of waves.

Joint models provide a flexible framework for characterizing the association between longitudinal and time-to-event processes and have been widely applied in biomedical research. However, fitting joint models can be computationally challenging for large-scale and complex biomedical data. This paper introduces the \proglang{R} package \pkg{FastJM}, which provides computationally efficient frequentist estimation for three classes of semiparametric joint models: joint models with a single longitudinal biomarker, joint models with multiple longitudinal biomarkers, and joint models with a single longitudinal biomarker with heterogeneous within-subject (WS) variability. Within an expectation–maximization framework, \pkg{FastJM} employs customized linear-scan algorithms to efficiently update the nonparametric baseline hazards, thereby addressing a major computational bottleneck in semiparametric joint modeling. The package also supports commonly used time-dependent latent association structures by integrating these algorithms with a landmark multivariate joint modeling framework. \pkg{FastJM} provides a unified interface for model specification, estimation, inference, visualization, dynamic prediction, and prediction performance assessment, including cross-validated time-dependent accuracy measures and time-independent concordance statistics. We describe the underlying methodology and software implementation and demonstrate the main functionality of \pkg{FastJM} through reproducible examples.

We propose a novel methodology for fitting decision trees to latent classes. The latent class analysis methodological literature has previously been focusing on logistic models of class membership given covariates, which has important drawbacks in the presence of complex interactions between covariates: logistic models are easily misspecified by omitting some interaction terms and complexity of interpretation increases rapidly with inclusion of higher-order interaction terms. Our proposed approach directly integrates decision tree modeling into the latent class analysis framework to model the covariate effects as an easily interpretable tree leading the eye through combinations of predictors to a final conditional classification. The novel methodology does not require any non-traditional assumptions for latent class models with covariates, is based on well-established routines of model estimation, and can be readily implemented in existing software. In the present paper, we focus on decision trees for binary covariates. We present the proposed approach, describe two tree pruning strategies, and provide a real-data illustration. Important extensions include generalizing the approach to multinomial and continuous covariates by developing more advanced within-predictor splitting procedures, and developing and evaluating alternative tree pruning strategies. The ultimate aim of this work is to initiate a novel research line in the latent class analysis methodological literature, and to facilitate the advancement of applied research via a novel data analytical tool.

Accurate small-area estimation of voter turnout for demographic subgroups is crucial for political analysis but methodologically challenging. Survey data suffer from over-reporting, non-representativeness, and non-ignorable non-response, while ecological inference (EI) from aggregate data is vulnerable to the ecological fallacy. We propose a Bayesian hierarchical framework that jointly integrates individual-level survey data, official aggregate turnout margins, and census cell counts. The framework comprises three models: a multilevel ecological inference (EI) model that extends MRP-style structure to aggregate data, a Poll-and-Margin (PM) model that jointly fits survey and margins in one coherent likelihood, and a Full Selection (FS) model that adds a Heckman-style selection correction under a random-contact assumption and uses informative priors for identification. In simulations with controlled selection mechanisms based on empirical census demographics, FS substantially improves over established benchmarks, particularly under strong selection bias or heavy censoring, but is sensitive to unmodeled selection noise. In an application to the 2023 Estonian parliamentary election, PM performs similarly to the current state of the art on available validation margins, with FS showing a modest improvement that depends on the informative priors.

Distributional mismatch between the data used to construct a statistical procedure and the population to which it is ultimately applied is pervasive in modern data analysis. We study covariate shift, a fundamental instance of this problem, and develop an adaptive importance-weighted model averaging method for prediction when labeled observations are available from a source distribution, whereas only unlabeled covariates are observed from the target distribution. Procedures fitted directly to the source sample generally optimize prediction risk under the source distribution and may therefore be suboptimal for target prediction. Importance weighting by the density ratio between the target and source covariate marginals provides a natural correction, but a small number of large density-ratio values can substantially inflate the variance of the resulting estimator in finite samples. We address this bias-variance trade-off by treating the degree of importance-weighting correction as a source of model uncertainty. Specifically, we construct a family of adaptive importance-weighted least-squares estimators by raising the estimated density ratio to a range of exponents, with the endpoints corresponding to ordinary least squares and standard importance-weighted least squares, and form a data-driven average over these candidates. Under model misspecification, the proposed model averaging estimator is shown to be asymptotically optimal relative to the infeasible best convex combination of the candidate estimators. Under correct specification, a diverging penalty is shown to make the selected weights concentrate near the ordinary least-squares endpoint. Simulations and a real-data application show that the proposed method achieves competitive target-prediction performance across the settings considered.

Bayesian statistics rests on a few familiar distinctions: frequentist vs. Bayesian, objective versus subjective probability, a model versus the data it is fitted to, a prior versus a posterior. Here, I use Donald Davidson's "third dogma of empiricism" to critique such distinctions in terms of scheme/content dualisms. With a single running example – estimating how often a football (soccer) team scores – I argue that a degree of belief is in itself objective, that model and data are better seen as parts of one belief system than as different kinds of things, and that Bayes' theorem describes relations within that system rather than data updating a model from outside. What matters instead is how belief systems are built, checked and rebuilt. To this end, I use Davidson's constructive programme of radical interpretation to suggest that a belief system should be not only coherent (as statisticians already require) but also charitable and triangulated, answerable to other people and to a shared world. I argue that these three constraints are exactly what good statistical practice – from eliciting priors to checking models – already tries to satisfy and talking about them explicitly could improve this practice.

The Weibull distribution is one of the most flexible continuous probability distributions used to model various failure rates and skewed data in reliability engineering, industry, weather studies and cancer studies. It is a common scenario in statistical inference that several Weibull populations share the same shape parameter, which also implies that they have the same coefficient of variation. While the inferential study of the common shape parameter is often considered for complete samples, the presence of censored data requires a separate investigation that has not received enough attention in the existing literature. Therefore, the focus of this article is on a comparative study of interval estimators for the common shape parameter and the common coefficient of variation under progressive type-II censoring using frequentist methods based on large-sample theory, variance estimates recovery, generalized pivots, and Bayesian inference. An optimal censoring scheme is also proposed to enhance the robustness of interval estimation. Numerical data analyses using a simulation study and a real carbon fiber strength data example are carried out for comparison, and the results recommend the intervals based on Bayesian and variance estimates recovery methods for their satisfactory performance.

Finding the optimal individualized treatment rule that maps individual characteristics or contextual information to treatment assignments has been extensively investigated in existing literature, with widespread practical applications. This paper considers the estimation of optimal treatment regimes within a semi-supervised data framework (exemplified by electronic medical record data). In such settings, only a tiny proportion of observations have observed outcome labels, owing to high labeling costs, time limitations, data privacy concerns, and other constraints, while covariates and treatment assignments are available for all study subjects. We develop a semi-parametric inference method for optimal treatment regimes, which leverages outcome- unlabeled samples with complete covariate and treatment information to enhance estimation efficiency. The proposed estimation framework consists of two key steps: first, flexible nonparametric imputation via single-index kernel smoothing; second, subsequent estimation of the optimal treatment regime based on concordance-assisted learning. We establish the consistency and asymptotic normality of our proposed estimators. Numerical simulation studies demonstrate that our method achieves higher efficiency and stronger robustness relative to fully supervised estimators under finite-sample settings. We further validate the practical value of our proposed framework using the MIMIC-III and ACTG175 datasets.

Forecast combining and forecast reconciliation for hierarchical and grouped time series have largely developed as separate research areas. This paper connects the two by developing a forecast combination framework for forecast reconciliation. For each bottom-level series, we construct a maximal linearly independent set of structured candidate forecasts from aggregation constraints, and show that combining these candidates and aggregating the resulting bottom-level forecasts is equivalent to standard unbiased linear reconciliation. Within this representation, we prove that mean-squared-error optimal combination weights exactly recover the widely used Minimum Trace (MinT) reconciliation. We further show that the optimal weight problem is separable across different bottom-level series, each yielding a Bates–Granger optimal forecast combination. This reveals MinT as a collection of optimal combinations over hierarchy-induced candidate forecasts. For finite-sample estimation, we propose a modular penalized framework that nests existing MinT variants and supports rich extensions including covariance shrinkage, weight penalization, and scalable series-wise separate estimation. Empirical results show that the framework is practically implementable, competitive with existing methods, and can improve accuracy while preserving coherence. Overall, the forecast combination perspective offers new interpretations of existing reconciliation approaches and provides a flexible basis for designing new methods.

High-dimensional functional panels consist of temporally ordered curves observed across many subjects and naturally exhibit heterogeneous structural changes. Under sparse subject-level break signals or opposite-signed shifts, traditional mean-aggregated CUSUM procedures may suffer noticeable power loss due to signal attenuation or cancellation induced by cross-sectional averaging. We propose a novel Energy–PE statistic, which combines subject-wise squared CUSUM energy aggregation with a generalized power-enhancement component. The energy aggregation preserves subject-level evidence under sign-heterogeneous changes, while the power-enhancement component improves sensitivity to sparse weak break signals. Under regularity conditions, we establish the asymptotic behavior of the proposed statistic. We further incorporate a latent group structure and an information-criterion-based clustering algorithm to estimate the unknown group number and membership for heterogeneous break points. Numerical studies and an intraday stock application demonstrate that Energy–PE controls size, improves power under sparse and sign-heterogeneous alternatives, and yields interpretable post-test summaries.

Receiver operating characteristic (ROC) and the area under the ROC curve (AUC) are widely used to evaluate the discriminative ability of biomarkers. In many clinical settings, however, diagnostic accuracy varies substantially across patient characteristics, and failure to account for such heterogeneity can lead to misleading conclusions. We propose a new semiparametric framework based on generalized additive models to estimate covariate-specific and covariate-adjusted AUC while allowing for nonlinear and interaction effects of covariates on biomarker performance. Our method accommodates both binary and multicategory disease status. We also establish the asymptotic properties of the proposed estimators. Simulations demonstrate favorable finite-sample performance. We illustrate the method using data from the Alzheimer's Disease Neuroimaging Initiative, where substantial heterogeneity in biomarker discrimination across covariates is observed.

Recovering geographic variation in survival requires separating spatial risk from patients' clinical characteristics, a problem complicated by prognostic covariates that are themselves spatially structured. We develop a nonparametric two-stage method for this separation. A clinical survival tree fit to the covariates alone supplies leaf Nelson-Aalen cumulative hazard residuals, transferring the censored survival structure to a clinically adjusted scale without imposing a functional form on the clinical hazard, and a second tree fit to these residuals on the coordinates recovers the spatial structure. Both stages are kernel dipole-splitting survival trees, so the resulting spatial risk map is piecewise constant, with sharp, possibly curved boundaries. We establish when the residuals recover the spatial signal: under a multiplicative frailty and an exogeneity condition, they are free of the clinical covariates given location and stochastically ordered by the frailty. An expansion of the frailty Laplace exponent quantifies their approximate exponentiality and identifies the leading remainder. These results assume consistency of the first stage rather than a model class, so the construction extends to other survival estimators. On the LeukSurv leukemia data the method agrees with a Bayesian Gaussian random field frailty about where risk is elevated while resolving sharp adjacencies the smooth surface averages away, and an unadjusted spatial analysis misattributes clinical variation to location. Simulations with known zones, including a sweep through graded violations of exogeneity, locate the point at which the two contributions cease to be separately identifiable, with the smooth benchmark degrading in parallel as that point is approached.

Despite growing interest in using LLMs to add robustness or reduce data-collection costs in survey experiments, their efficacy in conjoint design—an increasingly popular method in political science—remains underexplored. This paper addresses that gap by investigating whether synthetic agents can reproduce the multi-dimensional human preference patterns that conjoint is designed to capture. It replicates published conjoint studies and compares the results generated by synthetic agents with original human data along three dimensions: representational correspondence, inferential correspondence, and procedural stability. Our analysis evaluates the alignment of choice distributions as well as the statistical and substantive similarity of estimates, and the results are uneven across these dimensions and studies replicated. This implies that the validity of synthetic participants should be considered claim-dependent and hierarchical. Reproducing a figure or obtaining strong sign agreement is evidence of similar aggregate outputs, but not enough to support replacing human respondents. Our results suggest that the discipline as a whole must first map this innovation's boundaries across various levels before considering synthetic agents a robust substitute for human samples.

We study minimax-optimal designs and estimators for estimating the sample average treatment effect in finite population randomized experiments, where both design and estimator are unrestricted. For binary potential outcomes, we show this minimax risk is equivalent to the minimax risk $\rho_n^*$ of an estimation problem with $2$ unknown parameters. We leverage this reduction to establish a second-order risk expansion $\rho_n^* = n^{-1} - Cn^{-4/3} + o_n(n^{-4/3})$ for an explicit constant $C$ related to the Airy function. The minimax risk is attained by Bernoulli randomization with a nonlinear shrinkage estimator. Our results show that standard procedures such as complete randomization with difference in means are only minimax optimal up to first order in $n.$ We derive further results on admissibility of these procedures and discuss the practical implications of our results.

Machine learning (ML) has become an indispensable part of modern engineering design workflows. A crucial step in training an ML model is the selection of the loss function which can be systematically formulated via various techniques such as maximum likelihood estimation (MLE) and cross-validation . While MLE is one of the most popular, effective, and intuitive mechanisms for training ML models, it is brittle: if the assumptions underpinning it are not met, the trained ML model may generalize poorly. This brittleness affects even Gaussian processes (GPs) which are widely used in engineering design and are often (incorrectly) presumed to be very robust to overfitting. In this paper, we fundamentally evaluate the brittleness of MLE in the context of training GPs for probabilistic regression or classification tasks. We compare theoretically grounded metrics against MLE and propose practical solutions. Our extensive studies demonstrate the effectiveness of our solutions in downstream design tasks such as Bayesian optimization and provide a blueprint for practitioners to build accurate and robust GPs that can even outperform tabular foundation models in terms of prediction accuracy, uncertainty quantification, and inference cost. Our contributions are publicly available via GitHub at https://github.com/Bostanabad-Research-Group/GP-vs-TabPFN-vs-GPyTorch.

In the era of modern data science, the rapid proliferation of high-dimensional and functional datasets has fostered increasing interest in the investigation of paradigm shifts and structural breaks. Unlike classical univariate time series, structural changes in functional data need not occur simultaneously across the entire domain; instead, they may emerge locally, producing heterogeneous distortions across the underlying functional structure. The patterns of instability often exhibit sparsity, where it is not known a priori which specific parameters are undergoing a transition. However, in functional contexts, these shifts are often "localised". The difficulty lies in the high dimensionality of the parameter space, where the signal-to-noise ratio may be low for individual components, necessitating the aggregation of information across dimensions to detect a global change. This paper addresses the problem of detecting structural shifts in a functional time series from a Bayesian perspective. We have developed various novel methodologies that capture the inherent structural distortion in a sequence of random functions, both individually and simultaneously. The formulation of the problem is based on the state-space representation of a functional time series. Efficient Blocked Gibbs Sampling algorithms have been proposed to identify these locations accurately. Further, we demonstrate the effectiveness of our methods on several financial and temperature datasets.

Binary bipartite data often combine row and column heterogeneity with hierarchical interaction, as when students answer exercises organized by prerequisites or legislators vote on nested policy areas. Fixed-dimensional Euclidean volume grows only polynomially with radius, whereas branching hierarchies expand exponentially. Hyperbolic space matches this growth, and its rooted Gromov product measures shared ancestry. We propose the HypErbolic Latent Position model for bipartite Interaction data (HELPI), which places both object types in hyperbolic space while separating geometric interaction from additive propensities. With unrestricted main effects, a hyperbolic-distance predictor is likelihood-equivalent to a rooted Gromov-product predictor. Radial depth is absorbed by additive terms, leaving a root-invariant projected hierarchy signal as the identified geometric target. We establish identification and posterior contraction results and develop augmented variational procedures for binary outcomes. Simulations support recovery of hierarchical interaction and clarify weak-branch regimes. In Junyi Academy data, curriculum anchors yield an interpretable geometry with temporal prediction close to item-response benchmarks, while distinguishing attempted exercises from unconditional mastery. In U.S. House roll calls, flexible HELPI variants improve on additive and two-parameter logistic benchmarks and recover party and policy structure without using party labels during fitting.

Nonprobability samples (NPS) are attractive because they are less costly to collect, can provide substantially larger sample sizes, and may reach populations that traditional probability surveys do not. As response rates for traditional surveys fall, interest in NPS has grown rapidly within the field of survey statistics. These methods are especially relevant for small area estimation (SAE), where there is ever-present demand for estimates at fine geographic scales and detailed demographic domains. Despite rapid methodological development, there remains limited understanding of which approaches perform best under different conditions. In this paper, we review recent developments in NPS methodology, including the concept of data defect correlation (DDC) as a measure of data quality and as a tool for categorizing the various NPS methods. We then present a comprehensive simulation study that evaluates a range of NPS approaches under varying levels of DDC and extend several existing methods to the SAE setting.

Local polynomial smoothing is commonly used in non-parametric regression, but local linear derivative estimation still has a bias of order $O(h^2)$. This paper proposes an iterative data sharpening method to reduce the bias of derivative estimates while retaining the simplicity of local linear fitting. The method is based on two expectation operators: $L_0$, acting on the regression function, and $L_1$, acting on the first-order derivative. By repeatedly applying the residual operator $R=I-L_0$, a series of sharpened derivative estimates can be constructed. After $l$ sharpening steps, the bias order can be reduced from $O(h^2)$ to $O(h^{2l+2})$. For the Gaussian kernel, all sharpening coefficients equal 1, giving a simple closed-form single-bandwidth expression. Simulation experiments on three smooth test functions show that this method can significantly reduce the estimation bias while revealing a bias-variance trade-off.

Area-level small area estimation (SAE) models, such as the Fay–Herriot (FH) model, borrow strength across domains through covariates and random effects, but they can struggle when the relationship between the covariates and the outcome is spatially heterogeneous, that is, when it changes across the spatial domain of interest. We propose a spatially-clustered FH (SC-FH) framework that simultaneously (i) estimates cluster-specific regression coefficients and random effects variances and (ii) generates spatially coherent partitions of the geographical domain. Estimation maximizes a penalized likelihood that augments the FH likelihood with a Potts-type spatial cohesion term over the areal adjacency graph, through an efficient strategy that alternates between sequential label updates and closed-form FH updates within clusters. In simulation experiments run on the real geography of the application, the method recovers the latent regimes almost exactly whenever they are separated in the covariate–response space and improves prediction accuracy over the standard FH benchmark, with the spatial penalty acting as a stabilizer of both classification and estimation. An empirical application to the average standard output of farms in the Po Valley (Northern Italy) identifies two spatially compact production regimes with significantly different cluster-wise coefficients, and shows that the clusterwise predictor improves on the direct estimates while avoiding the over-shrinkage of the pooled model.