Papers

One of the major limits of kernel ridge regression (KRR) is that storing and manipulating the kernel matrix K_n for n samples requires O(n^2) space, which rapidly becomes unfeasible for large n. Nystrom approximations reduce the space complexity to O(nm) by sampling m columns from K_n. Uniform sampling preserves KRR accuracy (up to epsilon) only when m is proportional to the maximum degree of freedom of K_n, which may require O(n) columns for datasets with high coherence. Sampling columns according to their ridge leverage scores (RLS) gives accurate Nystrom approximations with m proportional to the effective dimension, but computing exact RLS also requires O(n^2) space. (Calandriello et al. 2016) propose INK-Estimate, an algorithm that processes the dataset incrementally and updates RLS, effective dimension, and Nystrom approximations on-the-fly. Its space complexity scales with the effective dimension but introduces a dependency on the largest eigenvalue of K_n, which in the worst case is O(n). In this paper we introduce SQUEAK, a new algorithm that builds on INK-Estimate but uses unnormalized RLS. As a consequence, the algorithm is simpler, does not need to estimate the effective dimension for normalization, and achieves a space complexity that is only a constant factor worse than exact RLS sampling.

Rejection sampling is a technique for sampling from difficult distributions. However, its use is limited due to a high rejection rate. Common adaptive rejection sampling methods either work only for very specific distributions or without performance guarantees. In this paper, we present pliable rejection sampling (PRS), a new approach to rejection sampling, where we learn the sampling proposal using a kernel estimator. Since our method builds on rejection sampling, the samples obtained are with high probability i.i.d. and distributed according to f. Moreover, PRS comes with a guarantee on the number of accepted samples.

We introduce and analyze a statistical estimator for Monge transport maps: solutions to the quadratic optimal transport problem in the Euclidean space. For absolutely continuous source measures, this map is uniquely defined as the gradient of a convex function, a result known as Brenier's theorem. Without absolute continuity, the problem is relaxed, maps are replaced by coupling measures, and optimal couplings are supported on the subdifferential of a convex function, called a Brenier potential. This potential is the basis for our transport map statistical estimator, for measures known only through finite samples. We construct an explicit optimal Brenier potential for the empirical problem, with a simple closed form expression based on dual solutions. Our transport map estimator is then defined as a suitable subdifferential selection, and can be extended to $\mathbb{R}^d$. We exhibit convergence rates for this estimator based on a new error bound for the quadratic optimal transport problem. In the semi-discrete setting, where the target measure is finitely supported, our estimator enjoys sharper convergence rates. Our methodology does not rely on smoothness or continuity of the target Monge transport map and requires no computation beyond primal-dual solutions of the empirical finite dimensional problem. Finally, using similar proof techniques, we provide novel convergence rates for empirical couplings.

Detecting copy number alterations (CNAs) from next-generation sequencing data remains challenging, particularly for short segments under noisy conditions. Existing segmentation methods often suffer from high false positive rates or fail to reliably detect short aberrations, especially in low-coverage data. In this study, we propose a modified tail-greedy unbalanced Haar (TGUHm) method that introduces a dual-thresholding strategy to improve segmentation accuracy. The proposed approach effectively suppresses spurious spikes while preserving sensitivity to both short and long CNA segments. Extensive simulation studies under Gaussian and heavy-tailed noise demonstrate that TGUHm consistently achieves higher true positive rates and lower false positive rates compared to state-of-the-art methods, including CBS, HaarSeg, and FDRSeg. In particular, the proposed method improves detection accuracy for short segments while maintaining competitive overall performance. Application to real cancer genomic data further confirms the practical utility of the method, revealing biologically meaningful CNAs associated with known cancer-related genes. These results suggest that TGUHm provides a robust and effective framework for CNA detection in challenging sequencing settings.

Classical ReLU-based Input Convex Neural Networks (ICNNs) are equivalent to the optimal value functions of Linear Programming (LP). This intrinsic structural equivalence restricts their representational capacity to piecewise-linear polyhedral functions. To overcome this representational bottleneck, we propose the SOC-ICNN, an architecture that generalizes the underlying optimization class from LP to Second-Order Cone Programming (SOCP). By explicitly injecting positive semi-definite curvature and Euclidean norm-based conic primitives, our formulation introduces native smooth curvature into the representation while preserving a rigorous optimization-theoretic interpretation. We formally prove that SOC-ICNNs strictly expand the representational space of ReLU-ICNNs without increasing the asymptotic order of forward-pass complexity. Extensive experiments demonstrate that SOC-ICNN substantially improves function approximation, while delivering competitive downstream decision quality. The code is available at https://anonymous.4open.science/r/SOC-ICNN-4B18/.

A nonparametric model using a sequence of Bernstein polynomials is constructed to approximate arbitrary isotropic covariance functions valid in $\mathbb{R}^\infty$ and related approximation properties are investigated using the popular $L_{\infty}$ norm and $L_2$ norms. A computationally efficient sieve maximum likelihood (sML) estimation is then developed to nonparametrically estimate the unknown isotropic covaraince function valid in $\mathbb{R}^\infty$. Consistency of the proposed sieve ML estimator is established under increasing domain regime. The proposed methodology is compared numerically with couple of existing nonparametric as well as with commonly used parametric methods. Numerical results based on simulated data show that our approach outperforms the parametric methods in reducing bias due to model misspecification and also the nonparametric methods in terms of having significantly lower values of expected $L_{\infty}$ and $L_2$ norms. Application to precipitation data is illustrated to showcase a real case study. Additional technical details and numerical illustrations are also made available.

Enhancing seismic fragility and risk assessment of nuclear power plants relies on accurate prediction of reactor building responses to seismic hazards, which can be further improved through dynamic analysis of high-fidelity finite element (FE) models. However, FE models often exhibit non-negligible discrepancies from actual structures due to various sources of uncertainty, necessitating FE model updating with rigorous quantification of associated uncertainties. This paper presents a GPU-accelerated latent space–based Bayesian framework for FE model updating of building structures. In the proposed framework, high-dimensional structural response data (e.g., time histories or frequency response functions) are projected into a low-dimensional latent space using a multimodal variational autoencoder (MVAE), thereby enabling efficient and tractable likelihood evaluation without explicit modeling in the original observation space. Once trained, the surrogate enables amortized inference, allowing posterior sampling to be performed without additional simulator evaluations. We specifically employ a sequential Monte Carlo (SMC) sampler, whose population-based formulation allows parallel evaluation of the approximate likelihood on GPUs, resulting in computational efficiency and robustness against multimodal and complex posterior distributions. The proposed framework is validated through both numerical benchmarking and experimental data from a shaking table test of a reinforced concrete building structure. The results demonstrate that the method accurately estimates structural parameters with well-quantified uncertainties, while achieving fast and efficient inference through GPU-based parallelization, and enabling robust inference even in the presence of sparse observations that induce multimodal and highly complex posterior distributions.

This paper studies expected performance and practical feasibility of the most commonly used classes of source-level likelihood-ratio (LR) systems when applied to a trace-reference comparison problem. The paper compares performance of these classes of LR systems (used to update prior odds) to each other and to the use of prior odds only, using strictly proper scoring rules as performance measures. It also explores practical feasibility of the classes of LR systems. The present analysis allows for a ranking of these classes of LR systems: from specific-source feature-based to common-source anchored or non-anchored score-based. A trade-off between performance and practical feasibility is observed, meaning that the best performing class of LR systems is the hardest to realise in practice, while the least performing class is the easiest to realise in practice. The other classes of LR systems are in between the two extremes. The one positive exception is a common-source feature-based LR system, with good performance and relatively low experimental demands. The paper also argues against the claim that some classes of LR systems should not be used, by showing that all systems have merit (when updating prior odds) over just using the prior odds (i.e. not using the LR system).

Most clinical prediction studies are developed from retrospective cohorts and reported as if all patient information were observed at once. In practice, clinicians face a more consequential question: when is there already enough information to stop testing and act? A later stage can produce a better-looking model and still fail to justify the added delay, burden, or invasiveness of further workup. We formulate sequential clinical prediction as an optimal-stopping problem under staged information, and illustrate the framework across four retrospective clinical datasets. The preferred stopping stage differed substantially by setting: sometimes fuller information justified waiting, whereas in other cases early or intermediate action was preferable. The key object is the patient-specific conditional risk trajectory: forward martingale structure represents coherent risk updating across stages, while reverse-martingale ideas describe information loss when a richer predictor is replaced by a simpler score. The results demonstrate that the best-performing model is not always the best stage for clinical decision-making.

Galaxy surveys provide noisy, incomplete, selection-affected population snapshots rather than complete evolutionary histories. We formulate galaxy evolution as joint inference of effective physical states, their population intensity, and cross-epoch dynamics. An effective state is a coarse-grained description defined by predictive adequacy and distinguished from measured variables and learned representations. The intrinsic population is a finite non-negative measure evolving through one-galaxy transport, source and removal terms, and nonlocal two-to-one merger jumps. A sub-probability observational kernel maps this population to survey catalogs and supplies the likelihood for state inference. Projected luminosity and stellar-mass functions generally do not obey closed dynamics and cannot identify the underlying channels alone. In a restricted IllustrisTNG proof of concept, a finite-time non-merger kernel and reduced merger operator are calibrated on training root merger trees and evaluated on held-out lineages. Dynamical propagation improves the normalized later mass–sSFR distribution relative to a static baseline, while the merger channel captures number loss from disappearing secondary progenitors. Under noisy and censored mock observations, the propagated prior improves posterior-stacked reconstruction of the latent population relative to a frozen-static prior. The framework provides an observation-directed, simulation-assisted route to galaxy state-and-dynamics inference.

We study stochastic multi-armed bandits in which the objective is a statistical functional of the long-run reward distribution, rather than expected reward alone. Under mild continuity assumptions, we show that the infinite-horizon problem reduces to optimizing over stationary mixed policies: each weight vector \(w\) on the simplex induces a mixture law \(P^w\), and performance is measured by the concave utility \(U(w)=\mathfrak U(P^w)\). For differentiable statistical utilities, we use influence-function calculus to derive stochastic gradient estimators from bandit feedback. This leads to an entropic mirror-ascent algorithm on a truncated simplex, implemented through multiplicative-weights updates and plug-in estimates of the influence function. We establish regret bounds that separate the mirror-ascent optimization error from the bias caused by estimating the influence function. The framework is developed for general concave distributional utilities and illustrated through variance and Wasserstein objectives, with numerical experiments comparing exact and plug-in influence-function implementations.

Optimization over the Stiefel manifold $\mathrm{St}(p,d)$, the set of $p \times d$ column-orthonormal matrices, is fundamental in statistics, machine learning, and scientific computing, yet remains challenging in the presence of non-convex, non-smooth, or black-box objectives. Existing methods largely rely on either convex relaxations or gradient-based Riemannian optimization, limiting applicability in derivative-free and highly multimodal settings. We propose BOOOM (Black-box Optimization Over Orthonormal Manifolds), a general-purpose framework for loss-function-agnostic optimization on $\mathrm{St}(p,d)$. The key idea is a global Givens rotation-based parametrization that maps the manifold to an unconstrained Euclidean angle space while preserving feasibility exactly. Building on this representation, BOOOM employs a structured, parallelizable, derivative-free search based on Recursive Modified Pattern Search, enabling systematic exploration through plane-wise rotations without requiring gradient information and facilitating escape from poor local optima. We establish a unified theoretical framework showing equivalence between angle-space and manifold optimization, transfer of stationarity, and global convergence in probability under mild conditions. Empirical results across diverse problems, including heterogeneous quadratic optimization, low-rank and sparse matrix decomposition, independent component analysis, and orthogonal joint diagonalization, among other widely studied settings, demonstrate strong performance relative to state-of-the-art methods, particularly in non-smooth and highly multimodal regimes. We further illustrate its practical utility through a novel supervised PCA formulation applied to metabolomics data in colorectal cancer.

The linear hazard regression model developed by Aalen is becoming an increasingly popular alternative to the Cox multiplicative hazard regression model. There are no methods in the literature for selecting among different candidate models of this nonparametric type, however. In the present paper a focused information criterion is developed for this task. The criterion works for each specified covariate vector, by estimating the mean squared error for each candidate model's estimate of the associated cumulative hazard rate; the finally selected model is the one with lowest estimated mean squared error. Averaged versions of the criterion are also developed.

Objectives: Large language models (LLMs) are increasingly used for clinical text summarization, yet structured methods to assess associated patient safety risks remain limited. Failure Mode, Effects, and Criticality Analysis (FMECA) provides a proactive framework for systematic risk identification but has not been adapted to LLM-generated clinical content. This study aimed to develop and validate a novel FMECA framework for the prospective assessment of patient safety risks in LLM-generated clinical summaries. Materials and Methods: An interdisciplinary expert panel (n = 8) developed a taxonomy of failure modes through literature review and brainstorming. Standard FMECA dimensions (occurrence, severity, detectability) were adapted into 5-point ordinal scales. The framework was applied to 36 discharge summaries from four patients, generated by an open LLM (GPT-OSS 120B) using real-world clinical data from the Geneva University Hospitals. Reviewers independently annotated the summaries across two rounds. Inter-rater reliability was assessed at failure mode, severity and detectability score levels. Usability and content validity were evaluated using an adapted System Usability Scale and structured feedback. Results: The final framework comprised 14 failure modes organized into categories. Inter-rater agreement improved between rounds, reaching moderate-to-substantial agreement for failure mode identification and good agreement for severity and detectability scoring. Usability was rated as good (mean SUS: 79.2/100), with high evaluator confidence. Discussion and Conclusion: This study presents the first FMECA-based framework for systematic patient safety risk assessment of LLM-generated clinical summaries. The framework provides a structured and reproducible method for identifying clinically relevant risks caused by these summaries.

Purpose: To investigate whether longitudinal changes in timing and magnitude of PA are associated with physical function (PF) in older women. Methods: Women from OPACH study with accelerometry at baseline and WHISH study W1 and W2 were included. Minute-level PA counts were averaged and smoothed as diurnal PA curves. Consecutive-visit change was modeled within periods (baseline–W1 and W1–W2) as a Riemannian deformation from earlier to later curves, with two-dimensional initial momenta characterizing timing and magnitude shifts. Multivariate functional principal component analysis (MFPCA) summarized coupled timing-magnitude patterns, and principal component (PC) scores and deformation energy were derived for each participant-period. Linear mixed-effects models related these features to RAND-36 PF, adjusting for baseline PF and covariates. Results: Mean PA deformation in both periods showed downward shifts in PA magnitude and temporal redistribution after 10:00. Top 15 PCs explained at least 90% of variability in both periods. PC1 captured diurnal PA increase/decrease, explaining 22.4% of variability for baseline–W1 and 20.8% for W1–W2. Among participants with complete PF scores and baseline covariates (N=1,157), higher PC1 scores, reflecting relative increase/maintenance of PA across day, were positively associated with PF (P<0.0001). Deformation energy, a metric for overall diurnal pattern change between visits, showed a significant interaction with period for PF (P=0.003), with a larger positive association during W1–W2 than during baseline–W1. Conclusions: In older women, longitudinal changes in diurnal PA accumulation were associated with PF. Riemannian deformation analysis identified clinically interpretable markers of PA pattern change that may capture functional-aging information not represented by conventional PA summaries.

We propose a unified probabilistic framework for modeling high-dimensional count tensors with excess zeros. Such data arise naturally in a variety of applications, including single-cell Hi-C experiments, where observations are represented as a third-order tensor indexed by genomic locus pairs and cells. We develop a zero-inflated Poisson tensor model that captures the underlying latent structure through a low-rank tensor decomposition. For single-cell Hi-C data, the proposed framework further accommodates heterogeneous cell populations through latent cluster structure and exploits the ordered nature of genomic loci via smooth latent representations. We develop a likelihood-based estimation procedure together with a Bayes-optimal classifier for distinguishing structural zeros from technical zeros, enabling principled false-zero detection, imputation, and uncertainty quantification. Theoretically, we establish identifiability of the proposed model and consistency of the proposed estimators. Simulation studies and analyses of single-cell Hi-C data demonstrate improved performance in false-zero detection, latent structure recovery, and clustering.

Interaction effects are ubiquitous in applied statistical modelling, yet their meaningful interpretation remains challenging. The classic Johnson-Neyman (JN) technique (Johnson and Neyman 1936) addresses this challenge for two-way interactions by identifying the regions of a moderator's range over which a focal effect is and is not statistically significant. The int3ract package for R implements the JN technique and its three-way extension (the Johnson-Neyman-Krause, or JNK, technique) for both frequentist and Bayesian models. The function JNK_freq() auto-detects models fitted via lm()/glm(), RSiena's siena(), or lme4's lmer()/glmer(), but can also be applied to multiplicative interactions from (virtually) any model family by supplying a coefficient vector and covariance matrix directly. For Bayesian Stochastic Actor-Oriented Models (SAOMs) estimated with multiSiena, or any model producing posterior draws, the function JNK_bayes() produces conditional posterior distributions. For two-way interactions, classic shaded confidence-band plots are created that visually demarcate significant and non-significant regions along the moderator range; three-way interactions yield colour-gradient heatmaps with optional crosshatch overlays for non-significant regions. The package is designed to encourage richer, region-specific reporting of interaction effects in place of the conventional single-slope spotlight approach.

In many longitudinal studies, a large number of variables are measured repeatedly over time, with substantial missing data. Existing methods, such as probabilistic principal component analysis (PPCA), are ill-equipped to handle such incomplete, high-dimensional longitudinal data, as they fail to account for the nested sources of variation and temporal dependency inherent in repeated measures. We introduce hierarchical probabilistic principal component analysis (HPPCA), a two-level probabilistic factor model that explicitly separates between-subject variance from time-varying within-subject dynamics. The within-subject latent factors are modeled by a Gaussian process. We develop an EM algorithm to handle missing data and flexible covariance kernels, accelerated by computationally efficient initializers. Simulation studies demonstrated that HPPCA robustly recovers model parameters subspaces and substantially outperforms both standard PPCA and multivariate functional PCA in imputation accuracy, even under heavy missingness and model misspecification. An application to the long COVID symptoms in the Researching COVID to Enhance Recovery adult cohort revealed that HPPCA effectively captured the data's hierarchical structure and its learned features significantly improved the prediction of clinical outcomes and the recovery of masked clinical records compared to exisiting methods.

Experimental designs that are minimax in the presence of model misspecifications have been constructed so as to minimize the maximum, over classes of alternate response models, of the integrated mean squared error of the predicted values. The theory to date has focussed almost exclusively on Least Squares estimates. Here we extend this theory to designs tailored for M-estimation of parameters, thus obtaining additional robustness against outlying, or otherwise erroneous, responses. We show that, subject to a minor change in a tuning constant, designs optimal for Least Squares remain so asymptotically for M-estimation, under squared error loss. We argue that even this minor change should be ignored, and the tuning constant chosen in an ad hoc but sensible manner that does not depend on which M-estimate is being employed. A rather surprising additional result is that our designs and estimates, derived under an assumption of i.i.d. errors, are also robust, in a minimax sense, against broad classes of correlation structures. We also investigate designs under a more robust loss function, less influenced by outlying responses than squared error.

Network meta-analysis (NMA) combines direct and indirect comparisons across a treatment network, but exact contribution decompositions that reproduce NMA estimates are lacking, especially for multi-arm trials with correlated contrasts. We develop a contrast-space projection formulation of NMA that expresses the estimator as a linear mapping of observed pairwise contrasts onto the consistency-constrained contrast space. Building on this representation, we define direct and indirect evidence through a canonical within-study reduction that removes algebraic redundancy and yields a unique, invariant study-level decomposition. The resulting covariance-aware weights exactly reconstruct the NMA estimator and can be further resolved into indirect path-level components. Under fixed effects, the same projection also represents the generalized Cochran Q decomposition into within-design heterogeneity and between-design inconsistency. The framework yields diagnostic and graphical tools, including forest plots, tension plots, and path-based visualizations. Applications to empirical networks illustrate how the approach provides a reproducible and interpretable account of evidence contributions in NMA.

We study the minimax sample complexity of multicalibration in the batch setting. A learner observes $n$ i.i.d. samples from an unknown distribution and must output a (possibly randomized) predictor whose population multicalibration error, measured by Expected Calibration Error (ECE), is at most $\varepsilon$ with respect to a given family of groups. For every fixed $\kappa > 0$, in the regime $|G|\le \varepsilon^{-\kappa}$, we prove that $\widetilde{\Theta}(\varepsilon^{-3})$ samples are necessary and sufficient, up to polylogarithmic factors. The lower bound holds even for randomized predictors, and the upper bound is realized by a randomized predictor obtained via an online-to-batch reduction. This separates the sample complexity of multicalibration from that of marginal calibration, which scales as $\widetilde{\Theta}(\varepsilon^{-2})$, and shows that mean-ECE multicalibration is as difficult in the batch setting as it is in the online setting, in contrast to marginal calibration which is strictly more difficult in the online setting. In contrast we observe that for $\kappa = 0$, the sample complexity of multicalibration remains $\widetilde{\Theta}(\varepsilon^{-2})$ exhibiting a sharp threshold phenomenon. More generally, we establish matching upper and lower bounds, up to polylogarithmic factors, for a weighted $L_p$ multicalibration metric for all $1 \le p \le 2$, with optimal exponent $3/p$. We also extend the lower-bound template to a regular class of elicitable properties, and combine it with the online upper bounds of Hu et al. (2025) to obtain matching bounds for calibrating properties including expectiles and bounded-density quantiles.

Geographic context is often consider relevant to motor insurance risk, yet public actuarial datasets provide limited location identifiers, constraining how this information can be incorporated and evaluated in claim-frequency models. This study examines how geographic information from alternative data sources can be incorporated into actuarial models for Motor Third Party Liability (MTPL) claim prediction under such constraints. Using the BeMTPL97 dataset, we adopt a zone-level modeling framework and evaluate predictive performance on unseen postcodes. Geographic information is introduced through two channels: environmental indicators from OpenStreetMap and CORINE Land Cover, and orthoimagery released by the Belgian National Geographic Institute for academic use. We evaluate the predictive contribution of coordinates, environmental features, and image embeddings across three baseline models: generalized linear models (GLMs), regularized GLMs, and gradient-boosted trees, while raw imagery is modeled using convolutional neural networks. Our results show that augmenting actuarial variables with constructed geographic information improves accuracy. Across experiments, both linear and tree-based models benefit most from combining coordinates with environmental features extracted at 5 km scale, while smaller neighborhoods also improve baseline specifications. Generally, image embeddings do not improve performance when environmental features are available; however, when such features are absent, pretrained vision-transformer embeddings enhance accuracy and stability for regularized GLMs. Our results show that the predictive value of geographic information in zone-level MTPL frequency models depends less on model complexity than on how geography is represented, and illustrate that geographic context can be incorporated despite limited individual-level spatial information.

Empirical Bayes methods are widely used for large-scale inference, yet most classical approaches assume homoscedastic observations and focus primarily on posterior mean estimation. We develop a nonparametric empirical Bayes framework for the heteroscedastic normal means problem with unequal and unknown variances. Our first contribution is a generalized Tweedie-type identity that expresses the Bayes estimator entirely in terms of the joint marginal density of the observed statistics and its partial derivatives, extending the classical Tweedie's formula to settings with heterogeneous and unknown variances. Our second contribution is to introduce a moment-generating-function representation that enables recovery of the full posterior distribution within the f-modeling paradigm without specifying or estimating the prior distribution. The resulting method provides a unified framework for point estimation, uncertainty quantification, and hypothesis testing while accommodating arbitrary dependence between means and variances. Simulation studies and real-data analysis demonstrate that the proposed approach achieves accurate shrinkage estimation and reliable posterior inference in heterogeneous data environments.

How can we monitor, in real time, whether one uncertain prospect has any upside over another? To answer this question, we develop a novel family of sequential, anytime-valid tests for stochastic dominance (SD), a classical and popular notion for comparing entire distribution functions. The problem is distinct from that of testing mean dominance, and it is particularly useful when comparing distributions with similar means or with ordinal outcomes. We first derive powerful, nonparametric e-processes that quantify evidence against the null hypothesis that one prospect is stochastically dominated by another. For first-order SD, these e-processes are based on mixtures of growth-rate optimal e-variables, yielding a test of power one that retains validity under continuous monitoring. We then generalize the approach to sequential testing for higher-order SD and other integral stochastic orders. Empirically, we find that the tests are competitive in power with classical, non-anytime-valid SD tests. Our real-world application examines a controversial phenomenon in baseball analytics, known as the "third-time-through-the-order (3TTO) penalty," viewed as a monitoring problem. We close by sketching the complementary problem of testing whether a prospect has a definite upside, formalizing conditions under which we can derive a powerful anytime-valid test.

Bayesian Optimal Experimental Design (BOED) provides a rigorous framework for decision-making tasks in which data acquisition is often the critical bottleneck, especially in resource-constrained settings. Traditionally, BOED typically selects designs by maximizing expected information gain (EIG), commonly defined through the Kullback-Leibler (KL) divergence. However, classical evaluation of EIG often involves challenging nested expectations, and even advanced variational methods leave the underlying log-density-ratio objective unchanged. As a result, support mismatch, tail underestimation, and rare-event sensitivity remain intrinsic concerns for KL-based BOED. To address these fundamental bottlenecks, we introduce an IPM-based BOED framework that replaces density-based divergences with integral probability metrics (IPMs), including the Wasserstein distance, Maximum Mean Discrepancy, and Energy Distance, resulting in a highly flexible plug-and-play BOED framework. We establish theoretical guarantees showing that IPM-based utilities provide stronger geometry-aware stability under surrogate-model error and prior misspecification than classical EIG-based utilities. We also validate the proposed framework empirically, demonstrating that IPM-based designs yield highly concentrated credible sets. Furthermore, by extending the same sample-based BOED template in a plug-and-play manner to geometry-aware discrepancies beyond the IPM class, illustrated by a neural optimal transport estimator, we achieve accurate optimal designs in high-dimensional settings where conventional nested Monte Carlo estimators and advanced variational methods fail.

Functional data analysis is concerned with the analysis of infinite-dimensional data functions. Functional principal component analysis (FPCA) is a key method to obtain finite-dimensional summaries. Consistency of FPCA has been theoretically established for sufficiently regular data functions. However, empirical evidence shows that FPCA can become severely inconsistent when the underlying functions are too rough. This paper provides the first theoretical explanation for this phenomenon. We propose a model that explicitly captures the roughness of functional data and allows us to quantify the resulting bias of FPCA, depending on the functional roughness. The model undergoes a phase transition marking the point at which FPCA becomes entirely uninformative. Based on these probabilistic results, we discuss diagnostic tests for informative principal components. As an additional contribution, we derive results on spectral statistics that may serve as a foundation for goodness-of-fit tests for rough functional data. Mathematically, our approach combines recent advances in random matrix theory and generic chaining with tools from FDA. We illustrate the effects of roughness on FPCA using simulations, as well as climate and environmental datasets.

Standard diffusion models are flexible estimators of complex distributions, but they do not encode causal structures and therefore do not by themselves support causal analysis. We propose a causality-encoded diffusion framework that incorporates a known directed acyclic graph by training conditional diffusion models consistent with the graph factorisation. The resulting sampler approximately recovers the observational distribution and enables interventional sampling by fixing intervened variables while propagating effects through the graph during reverse diffusion. Building on this interventional simulator, we develop a resampling-based test for directed edges that generates null replicates under a candidate graph. We establish convergence guarantees for observational and interventional distribution estimation, with rates governed by the maximum local dimension rather than the ambient dimension, and prove asymptotic control of type I error for the edge test. Simulations show improved interventional distribution recovery relative to baselines, with near-nominal size and favourable power in inference. An application to flow cytometry data demonstrates practical utility of the proposed method in assessing disputed signalling linkages.

In both observational studies and randomized trials, post-treatment events such as dropout, nonadherence, and truncation by death occur frequently. In some studies, conditioning on post-treatment variables is a deliberate strategy to isolate particular treatment effects on the outcome. However, naive comparisons of outcomes conditional on post-treatment events generally lack a causal interpretation, even when treatment is randomly assigned. There exist causal estimands that account for post-treatment events, including survivor average causal effects and conditional separable effects, but identification usually requires measurement of common causes of the post-treatment event and the outcome. In this article, we show that such measurements are not always necessary. Conceptually, what we require is that the treatment and other unmeasured causes of the outcome generate the post-treatment event through "independent mechanisms". Then, conditional separable effects and survivor average causal effects can be identified without adjustment for common causes of the post-treatment event and the outcome. We illustrate the results in studies with truncating events, differential nonadherence, and the birth weight paradox.

Diffusion-based generative models have reformed generative AI, and also enabled new capabilities in the science domain, e.g., fast generation of 3D structures of molecules. In such tasks, there is often a symmetry in the system, identifying elements that can be converted by certain transformations as equivalent. Equivariant diffusion models guarantee a symmetric distribution, but miss the opportunity to make learning easier, while alignment-based simplification attempts fail to preserve the target distribution. In this work, we develop quotient-space diffusion models, a principled generative framework to fully handle and leverage symmetry. By viewing the intrinsic generation process on the quotient space, the exact construction that removes symmetry redundancy, the framework simplifies learning by allowing model output to have an arbitrary intra-equivalence-class movement, while generating the correct symmetric target distribution with guarantee. We instantiate the framework for molecular structure generation which follows $\mathrm{SE}(3)$ (rigid-body movement) symmetry. It improves the performance over equivariant diffusion models and outperforms alignment-based methods universally for small molecules and proteins, representing a new framework that surpasses previous symmetry treatments in generative models.

Mendelian randomization is a powerful tool for causal inference in observational studies. The two-sample summary-data design, which estimates genetic associations with exposures and outcomes in separate cohorts, is the most widely used Mendelian randomization approach in large-scale genomic studies. However, this approach relies on a strong assumption of population homogeneity across the two samples. In practice, available samples often differ in ancestry, demographics, socioeconomic factors, covariate adjustment, and measurement protocols. Violations of the homogeneity assumption can bias causal effect estimates and undermine the credibility of Mendelian randomization findings. We introduce a robust, model-free Mendelian randomization framework that directly addresses population heterogeneity in the two-sample summary-data setting. Our method avoids parametric assumptions about population differences and is designed to address real-world challenges, including measurement error, weak instruments, and pleiotropy. We show that the proposed estimator is consistent and asymptotically normal under heterogeneous designs, and may offer efficiency gains over the classic estimator even in homogeneous settings. Through numerical simulations and a real data analysis for estimating the causal effect of body mass index on high-density lipoprotein cholesterol across ancestrally diverse populations, we demonstrate the practical utility, stability, and robustness of our approach.