Papers

Given a bivariate copula $C$, the bivariate copula $\Psi(C)$ is defined as the Markov product of $C^T$ and $C$. In this paper we consider the diagonal section of $\Psi(C)$ under the additional assumption that its value at a single point is known. We determine the exact upper and lower bounds for the diagonal section of $\Psi(C)$, given that $\Psi(C)(a,a)=d \in [a^2,a]$, where $a \in [0,1]$ is a fixed point. Furthermore, we describe some additional properties of the transformation $\Psi$, such as its behavior under ordinal sums and the survival operation. We also study checkerboard copulas in the range of $\Psi$ and matrices associated to them.

Gaussian process solvers for linear partial differential equations (PDEs) return posteriors rather than point estimates, but they condition only on equality information: the PDE residual, boundary data, and measurements. Many solutions of interest also obey inequality constraints, including nonnegativity, monotonicity, and absence of arbitrage, which equality conditioning cannot enforce. We prove that conditioning a Gaussian process prior jointly on linear equality observations and a finite-dimensional inequality event is equivalent to conditioning sequentially in terms of Borel measures on the solution space. The physics-informed posterior is therefore a valid prior for the inequality update at virtual points. The combined law is non-Gaussian yet computable through truncated-Gaussian moments, and neither the PDE conditioning nor the inequality update requires modification. We evaluated the construction on a Poisson benchmark, option-price surfaces, and air-quality fields. Under incomplete boundary data, a monotonicity constraint lowered the error of a two-dimensional Poisson reconstruction by $63\%$. On a real S&P 500 option chain thinned to $20$ retained quotes, no-arbitrage constraints reduced held-out extrapolation error by $22\%$ and cut the calendar-arbitrage rate from $10\%$ to $1\%$. On four real air-quality snapshots, nonnegativity reduced a Gaussian moment-based sub-zero-tail proxy by $74\%$ at the sparsest monitor count. Changes in mean accuracy and constraint diagnostics differed across data regimes. The probability of the noisy virtual constraint event under the equality posterior measures agreement with the specified bounds before the update. The construction turns admissibility from a post-processing step into conditioning information for probabilistic PDE solvers.

What can an adversary recover from a model trained with DP-SGD? Its $(\varepsilon,\delta)$ guarantee does not answer this question directly. We derive a distortion bound that can be computed before training and applies even to an adversary who observes the parameter trajectory during training, knows every other training record and the population distribution from which the training records are drawn. The analyst may choose the representation in which error is measured, allowing the bound to target the aspects of reconstruction that are privacy-sensitive in a given application. Our information-theoretic framework separates two contributions: the information about a record revealed by training, measured through conditional mutual information, and the uncertainty already present in the data distribution, measured through source entropy. Rate-distortion theory combines them into a lower bound on the expected error of every reconstruction attack. For DP-SGD, the leakage bound depends only on the sampling rates and noise multipliers, while the source entropy is lower-bounded from a sample with finite-sample confidence. An orthogonal block construction keeps the resulting bound informative for high-dimensional images. We instantiate the framework in pixel space and in LPIPS feature space, obtaining, to our knowledge, the first lower bound on reconstruction distortion measured directly by LPIPS. On a range of image datasets the bounds establish non-trivial unavoidable expected reconstruction error within heuristically recommended privacy-budget ranges and beyond. Experiments compare the bounds with strong reconstruction attacks and identify the principal sources of the remaining slack.

Prior-data fitted networks (PFNs) are trained to approximate the Bayesian posterior predictive distribution (PPD) of a synthetic prior, yet results about the approximation quality are surprisingly thin. Taking an empirical risk minimization perspective, we decompose the deployment error of PFNs into prior, context, and amortization error, and measure the latter against a reference PPD under a realistic prior. Results show that PFNs generally recover the targeted PPD very well. Nevertheless, we find that parts of the residual amortization error persist, especially in rare settings under the prior pretraining distribution, even when scaling the model size, which manifests in disagreement between independently pretrained PFNs and leaves uncertainty about the predictor at inference time. As a remedy, we study multiple types of Bayesian model averages and find that these can consistently reduce this uncertainty and improve predictions up to full-size TabICLv2 models.

We present a simulation-based method for the estimation, sensitivity analysis, and optimization of distortion risk measures (DRMs), a class of performance measures that capture a wide range of risk attitudes through an underlying distortion function. Existing simulation approaches often approximate DRMs using numerical partitions or empirical distributions. We show that, when the distortion function is piecewise polynomial, the nonlinear structure of a DRM can be handled exactly, without introducing distortion discretization error, using only simulated system outputs. This representation enables the development of a computational method that integrates finite-difference techniques for sensitivity analysis and for solving a broad class of differentiable risk-based optimization problems. Under suitable conditions, we establish strong consistency and convergence rates for DRM estimation, sensitivity analysis, and optimization. We further extend the approach to non-polynomial distortion functions through polynomial approximation. Numerical experiments on synthetic and calibrated financial examples demonstrate the effectiveness and robustness of the proposed approach across a range of applications.

Sensitivity analysis asks how strong unmeasured confounding would have to be before an observational conclusion changes. Established methods answer this question using different quantities, including partial $R^2$ values, correlation products, and risk ratios. In this article, we present confoundvis, an R package that organizes these results through a common object, plotting functions, and reports. It calculates partial $R^2$ and impact-threshold paths from fitted linear models, imports results from sensemakr, konfound, and EValue, and provides measured-covariate benchmarks for the OLS frameworks. We demonstrate the workflow with 5,735 patients from the right heart catheterization study. The adjusted linear probability coefficient is 0.0634, the point-null robustness value is 0.0608, and the conditional significance impact threshold is 0.0376. The largest adverse leave-one-term-out measured-covariate product is 0.0090. Because the benchmark and threshold use different conditioning sets, their ratio is descriptive. A separate modified Poisson model gives a risk ratio of 1.103 and an E-value of 1.439. The contribution is a software workflow that keeps each method's units and decision rule clear. The package is available from CRAN.

Neural networks (NNs), while highly flexible and versatile, are often poorly calibrated. Post-hoc calibration methods adjust the outputs of an already trained model using a separate calibration dataset and additional calibration parameters. This process adds overhead that requires separate optimization of the calibrator and subsequent model selection. In this paper, we show that standard (post-hoc) calibrators can be trained end-to-end with the backbone via a bilevel Stackelberg formulation. This formulation eliminates the need for a separate calibration set or a dedicated calibration loss by jointly optimizing accuracy and calibration using cross-entropy, introducing only one additional hyperparameter: the calibrator learning rate. Although our framework supports calibrators ranging from temperature and matrix scaling to NN layers, our results show that using the final classification layer itself as the calibrator yields accurate, well-calibrated models with no additional parameters. Our approach outperforms SOTA calibrators with minimal added complexity, while also avoiding the multiplicative overhead of separate calibration components in settings such as multi-task learning.

Linearizing softmax attention offers a route to efficient Transformer inference, but preserving pretrained attention remains challenging. We identify an important source of approximation error overlooked by standard linearization methods: a small number of outlier dimensions contribute disproportionately to query-key dot products and account for a large fraction of the approximation error. Motivated by this heterogeneity, we propose Outlier-Aware Attention Linearization (OAL), a dimension-aware method for constructing quadratic attention approximations. Starting from a second-order Taylor approximation, we introduce a Hadamard-product formulation that makes individual query-key contributions explicit. This formulation relaxes globally shared coefficients to support dimension-specific linear terms and cross-dimensional quadratic interactions while preserving separable evaluation. To retain this flexibility with fewer parameters, OAL groups dimensions with similar contribution scales to share coefficients, while allowing distinct treatment of the few outliers. The resulting positive quadratic kernel is calibrated against pretrained attention outputs and admits linear-time evaluation with respect to sequence length for fixed head dimensions. On Qwen2.5-0.5B, OAL reduces frozen-attention output relative MSE by approximately 8.9% compared with DiJiang. Replacing middle-layer attention modules and jointly fine-tuning the kernel coefficients with LoRA further lowers WikiText-103 perplexity from 31.71 to 27.50 compared with a learned global quadratic, while improving accuracy on PIQA, HellaSwag, and ARC-E. These results demonstrate the value of outlier-aware dimension grouping for both attention approximation and pretrained-model conversion. The code is available at https://github.com/wuqianyue09/OAL

Latent field models are essential for spatial statistics, but the posterior distribution of the latent field, even with a Gaussian process (GP) prior, is analytically intractable whenever the likelihood is non-Gaussian. Existing approximations are either too slow or too restrictive for large non-Gaussian datasets. We present a scalable variational inference framework that mirrors the sparse nearest-neighbor autoregressive structure of the prior in the variational family, uses normalizing-flow conditionals for expressiveness, and amortizes all variational parameters via a single shared neural conditioner. Doubly stochastic optimization with reduced ancestor sets achieves polylogarithmic complexity per iteration.

Reinforcement learning with verifiable rewards (RLVR) has substantially improved LLM reasoning, yet growing evidence suggests it mainly performs distribution sharpening: concentrating probability mass on reasoning paths the base model already produces. Recent work therefore seeks RLVR's gains without needing reward signals, by sampling from the power distribution $\pi_\alpha \propto \pi_{\text{base}}^\alpha$ via Markov chain Monte Carlo or sequential Monte Carlo. We show that these samplers fail to converge to $\pi_\alpha$ under realistic compute budgets, and that their gains are largely recreated via simple Best-of-$N$ with appropriate grading . This suggests they succeed by selecting high-probability sequences, not by approximating $\pi_\alpha$, raising a natural question: is the power distribution itself a worthwhile target? We test this hypothesis directly by fine-tuning the base model to match $\pi_\alpha$, sidestepping sampler convergence issues entirely. Despite approximating the target far more faithfully than existing samplers, this strategy yields gains that are noisy and specific to certain models and sampling tasks. Together, our results suggest that the value of power sampling lies not in faithfully targeting $\pi_\alpha$, but in the simpler goal of extracting sequences the base model already deems likely.

Circular data recorded over time, such as wind direction, are increasingly observed in intervals rather than as single point-valued observations, yet no existing framework accommodates data that are simultaneously circular, interval-valued, and temporally dependent. Time series models for point-valued circular data, including recently proposed score-driven approaches, assume that each observation is a single angle. This paper addresses this gap by proposing a dynamic conditional score framework for interval-valued circular time series. The lower and upper boundary sequences of the observed intervals are treated as circular time series governed by score-driven dynamics, which preserves the periodicity of the circular domain and the dependency of boundaries, in which the boundaries are coupled through a toroidal distribution and a jointly derived score. Maximum likelihood estimation is discussed, and their finite-sample performance is evaluated through simulation. An application to interval-valued wind direction data illustrates the practical value of retaining the interval structure within a fully rotationally invariant circular framework.

The external validity of a study is the extent to which its findings generalize to unseen experimental conditions, such as datasets. We formalize the results of a study together with its research question, and define external validity as the probability that two independent and same-sized replications of the study produce results whose maximum mean discrepancy (MMD) is below a given tolerance. This approach, centered on the choice of a space of results and of an appropriate reproducing kernel for the MMD, is flexible enough to be deployed on multiple different types of results and research questions. Additionally, we establish an asymptotic power-law relation between the size of a study and the quantiles of the MMD, which we use to estimate the desired size of an externally valid study from a small-scale preliminary study. We deploy our framework to rankings with ties and discuss kernels suited to different research questions. Experiments on existing case studies and on a simulated newly designed study illustrate the applicability of the framework. We release our code as part of the genexpy Python module.

We discuss three simple connections between generative diffusion and optimal transport. First, the joint law of the data and of their Ornstein-Uhlenbeck noising is, at every time, the solution of an entropic optimal-transport problem, at a temperature set by the noise level. Backward diffusion can consequently be read as an annealing which gradually decreases the entropic regularization. Second, we study the speciation transition–how the modes of the data distribution emerge during generation–and compare how it happens in diffusion and flow matching with respect to optimal transport. Using Gaussian mixtures as high-dimensional example, we show that the dynamical instability through which diffusion and flow matching separate the modes becomes, for optimal transport, an exponentially thin boundary layer at the beginning of generation, which is rounded by entropic regularization. Third, considering a generic forward dynamics with a general confining potential, we show that the backward process driven by the score also satisfies a drift-based variational principle whose minimal cost is the entropy difference between initial and final distributions.

Correlation matrices combine nonlinear positive-definiteness constraints with a simultaneous row-column permutation symmetry. Their generation therefore involves two distinct problems. First, every output must remain on the correlation elliptope, irrespective of how its variables are ordered. Second, the generator must learn the domain-specific distribution occupied by empirical data inside that valid space. Existing financial image generators mix these problems: they learn ordered matrix entries and repair invalid draws afterward. We introduce CorrFlow, which separates them. A generalized Fisher transform maps full-rank correlation matrices bijectively to unconstrained coordinates, guaranteeing valid outputs without projection. A graph flow-matching network learns the empirical distribution in those coordinates and is permutation equivariant: reordering the variables only reorders the output. We compare a conditional CorrGAN reimplementation, image diffusion, and CorrFlow on $80\times80$ Taiwan equity matrices using five seeds per neural method. CorrFlow produces 100% valid draws without repair and obtains feature sliced-Wasserstein distance $.712\pm.141$, versus $4.339\pm.252$ for CorrGAN and $4.610\pm.191$ for image diffusion after their invalid outputs are projected back to the correlation space. It also improves this distance over historical lookup, while lookup retains higher local coverage. These results support the combined construction under the evaluated training budget; they do not establish superiority over fully tuned image models or isolate the contribution of each component.

Thyroid nodules affect a large proportion of the general population. Although most are benign, a non-negligible fraction is malignant. Ultrasound imaging is the primary diagnostic modality, and accurate nodule segmentation is a prerequisite for CAD systems to classify malignancy and assess metastatic risk.\\ In this work, we propose a thyroid nodule segmentation framework that integrates a DeepLabV3+, channel-attention modules and post-hoc conformal prediction, providing both accurate mask and pixel-level uncertainty estimates. We address thyroid ultrasound nodule segmentation as the first stage of a pipeline that also includes malignancy classification and metastatic risk prediction. We adopt DeepLabV3+, whose Atrous Spatial Pyramid Pooling (ASPP) captures multi-scale context, with an EfficientNet-B5 backbone for its accuracy-efficiency trade-off. Additionally, we test three channel-attention modules: Squeeze-and-Excitation, CBAM and Grouped Channel Attention (GCA), paired with two loss formulations (Dice-BCE and Boundary).\\ Since binary masks do not indicate prediction confidence, we wrap the trained model with a split conformal prediction layer that produces pixel-level uncertainty maps with statistical coverage guarantees.\\ We evaluate all configurations on the TN3K public benchmark using a fixed split and directly compare them with state-of-the-art methods. The plain backbone already improves upon previously reported results, with CBAM + Dice-BCE emerging as the best overall configuration.

We analyze distributional drift in nonstationary weighted risk minimization by decomposing the excess risk into estimation, drift, and approximation components. The drift error is commonly bounded through discrepancy-type quantities that control changes of the risk over classes of predictors. We instead relate the drift error directly to the movement of the best predictor within the model class and analyze its behavior in parametrized regression models. Our analysis shows that discrepancy-type quantities can be overly conservative: unlike our oracle drift, they may scale at a larger order than the resulting excess-risk contribution and can respond to changes in the risk surface that do not affect the optimal predictor. For linear regression, we characterize this gap sharply and combine the resulting drift bounds with estimation guarantees. We further extend the analysis to nonparametric regression.

Uncertainty quantification for electric vehicle (EV) charging demand is essential to turn the increasing pressure of EV adoption on power systems into a source of flexibility. We study a large-scale dataset of public charging stations in Scotland, where demand exhibits a hierarchical structure, organized across multiple aggregation levels, from individual charging points to regional and national scales. The charging infrastructure is continuously evolving, with new stations deployed and existing ones occasionally removed, which calls for forecasting and uncertainty quantification methods tailored to this evolving structure. To address these challenges, we combine hierarchical forecasting with recent advances in conformal prediction with missing data. The conformal procedures considered produce valid prediction intervals at every node of the hierarchy and the algorithm CP-MNR-Nested successfully accounts for the time-varying hierarchical structure of the data.

Interferometric Synthetic Aperture Radar (InSAR) offers a low-cost, wide-area solution for near continuous infrastructure monitoring. However, translating InSAR observations into automated anomaly warnings remains challenging because displacements are relative, affected by uncertainty, and acquired along line of sight (LOS). We introduced a probabilistic warning framework that accounts for these factors. LOS displacements were modelled through a Gaussian distribution. A limit state function was defined to describe the displacement exceedance, and the time-dependent probability of exceedance was evaluated. Instead of a fixed threshold, warnings were identified through abrupt increases in exceedance probability. The framework was validated using three case studies of quay walls in Amsterdam, providing timely anomaly warnings and issued no warnings for the non-collapse case. Results demonstrate the potential of a probabilistic interpretation of InSAR time series for deformation anomaly warning.

When a language model is distilled from a teacher that shares its base parameters, the student can inherit behavioural traits of the teacher—including misalignment—even from data that is semantically unrelated to the trait (Cloud et al., 2026). All existing defences require knowing the trait. We remove this requirement and prove that a subliminally transferred trait can be recovered from the parameter updates alone, without supervision, by exploiting a single asymmetry: the trait is expressed consistently across diverse datasets generated by the teacher, whereas task-specific and shared capabilities vary. The resulting estimator returns a weight-space direction that can be amplified to discover and interpret the trait or projected out to suppress it. Applied to the subliminal-learning benchmarks of (Cloud et al., 2026), our method recovers the transmitted animal-preference traits without any trait labels and suppresses their expression with almost no impact to task performance.

This paper considers the problem of estimating the ratio of variance components of a balanced one-way random effects model. Existing estimators of ratio of variance, such as the ANOVA and maximum likelihood estimators, can take negative values and are inadmissible under quadratic loss. We therefore study the problem under a general class of scale-invariant loss functions. We characterize the class of scale-equivariant estimators of ratio of variance and obtain the best equivariant estimator. We then derive a class of estimators that dominate the BEE and extend this class to obtain further improved estimators. We propose a smooth improved estimator of ratio of variance. For interval estimation of ratio of variance, we derive an asymptotic confidence interval, bootstrap confidence intervals, a generalized confidence interval and a highest posterior density credible interval. The performance of the proposed improved point estimators and confidence interval procedures is assessed through an extensive simulation study. Finally, we provide a real data analysis for implementation purposes.

In spatial data analysis, global measures of spatial correlation provide useful insights into the general spatial patterns, but fail to capture local variations. Local Indicators of Spatial Association (LISA) address this limitation by identifying local clusters and spatial outliers, offering valuable applications in many fields of research such as public health, urban planning and environmental monitoring, to name only a few. However, the application of LISA tools in hypothesis testing is undermined by the presence of the multiple testing problem which dramatically increases the risk of committing Type I errors especially when dealing with very large datasets. Standard corrections procedure (like, e. g., Bonferroni and Benjamini-Hochberg procedures) fail to consider the spatial dependence among the various LISA tests and, furthermore, they are overly conservative. This paper proposes an adjustment that, by incorporating the notion of the "effective sample size", improves statistical inference balancing false discovery rate control and statistical power.

The BROJA partial information decomposition provides an operationally motivated way to separate the information that two sources carry about a target into unique, redundant, and synergistic components. For continuous and high-dimensional variables, however, computing the BROJA decomposition requires solving a difficult optimization problem. Previous work showed that, when each source-target pair is jointly Gaussian, the problem can be reduced to finite-dimensional numerical optimization. Here, we derive a closed-form solution to this optimization problem, enabling direct computation of the BROJA decomposition in the general multivariate Gaussian setting. The resulting closed-form expression enables speedups exceeding two orders of magnitude for high-dimensional Gaussian BROJA PID. Code is available at https://github.com/brainsandmachines/Eigen_PID.

This paper deals with the problems of decision theoretic estimation of the scale parameter, $\sigma$, and the entropy (Shannon or Renyi) in general scale mixture models with unknown location. Employing Kubokawa's approach, under convexity of the loss $L(t)$ either in $t$ or in $lnt$, improved estimators of $\sigma$ are derived. If $L(t)$ satisfies an additional condition, the problem reduces to estimating $\sigma$ in the degenerate (scale mixture) model. As the entropy is a linear translation of $ln\sigma$, its estimation under a loss $L_0(u)$ amounts to estimating $\sigma$ under $L(t)=L_0(lnt)$, with the quadratic loss $L_0(u)=u^2$ and Brown's loss $L(t)=(lnt)^2$ being a typical case. Thus, in the case that $L_0(u)$ is convex improved estimators of the entropy are directly obtained from respective estimators of $\sigma$. The results are applied to scale mixtures of normal and exponential distributions, including the multivariate $t$ and multivariate Lomax distributions.

Quantitative genetics is the study of continuous traits and the statistical analysis of phenotypic variation, that is, the observable differences in physical traits among individuals within a population. If family relationships between individuals are known, phenotypic variation can be partitioned into genetic and environmental components. Some traits, such as growth trajectories, can be described as a function of a continuous index and are known as function-valued traits. However, function-valued traits from different individuals often vary not only in their shape (amplitude) but also in the speed at which the biological process develops (phase). In this work, we present and compare two approaches to account for amplitude and phase variations in estimating the additive genetic and environmental components. The first one follows from a conventional framework where warping functions remove phase variations, and then the aligned curves are analysed. The second approach adopts the perspective of looking at the aligned curve and the warping function as a joint functional object. In both cases, we propose a functional mixed-effects framework that uses functional principal components as the data-driven basis for genetically correlated data. The proposed procedure is applied to simulations and a real-case study of Tribolium castaneum growth trajectories.

A realized tail-risk conditional autoregressive expectile (RTCARE) model is proposed for joint forecasting of value-at-risk (VaR) and expected shortfall (ES). Realized tail measures such as the realized quantile and the realized ES, constructed from intraday high-frequency returns, are incorporated into the dynamics and measurement equations of the conditional expectile, allowing information about lower-tail risk to enter the forecasts. The parameters are estimated by Markov chain Monte Carlo, and we evaluate the resulting VaR-ES forecasts in numerical examples and in an empirical application using U.S. stock-return data. Across the considered tail levels, the RTCARE models deliver competitive and often superior forecasting performance relative to existing alternatives, with stable results across Monte Carlo replications and rolling forecast evaluations.

Compression schemes provide a general framework for deriving data-dependent risk bounds for learning algorithms. For learning rules based on preferent compression, existing results upper-bound the risk in terms of a sample-dependent complexity $\boldsymbol{k}$, defined as the cardinality of the compression augmented by the training examples on which the reconstructed hypothesis is inappropriate. These results, however, yield a margin of order $O(\sqrt{\ln N/N})$ above $\boldsymbol{k}/N$. In this paper, we establish that the logarithmic factor is not intrinsic. We prove that, for every fixed confidence level $1-\delta$, the risk is bounded with probability at least $1-\delta$ by $$ \frac{\boldsymbol{k}}{N} + O\left(\frac{1}{\sqrt{N}}\right), $$ uniformly over all possible values of $\boldsymbol{k}$. Thus, risk estimates for learning algorithms based on preferent compression attain the same law-of-large-numbers rate as the empirical estimate of the probability of a fixed event. We also show that the leading term $\boldsymbol{k}/N$ is necessary and investigate its relation to the empirical risk. In particular, we identify regimes in which the empirical risk retains the $O(1/\sqrt{N})$ convergence rate, as well as regimes in which it converges more slowly or ceases to be consistent, while the complexity-based estimator preserves the $O(1/\sqrt{N})$ rate.

Time series causal discovery — the process of detecting causation between time series — has its roots in signal processing research from over a half century ago. It has since received active and widespread interest, evolving into different forms across domains such as statistics, dynamical systems and machine learning. Entering this fragmented and constantly evolving research landscape can be daunting. This tutorial provides a unified overview of time series causal discovery across research communities. We introduce the key principles invoked by modern discovery methods, and we discuss their limitations and assumptions. We also discuss how discovery methods are evaluated, benchmarked and fine-tuned. We further discuss how causal discovery for time series is becoming crucial for downstream tasks in agriculture and Earth science, in healthcare, and in artificial intelligence systems. We believe that this tutorial will serve as a valuable resource in navigating the deluge of recent causality research, as well as an invitation to signal processing audiences to contribute to causal problems. To further help readers, we host a companion repository to link to useful code packages and resources at https://github.com/CHAI-UK/Awesome-Time-Series-Causal-Discovery.

We develop a model-free, noise-aware Bayesian clustering framework that constructs a martingale posterior over clusterings from a predictive Bayesian perspective. This avoids explicit likelihood-prior model specification, expensive MCMC sampling, and mode-seeking variational approximations, while accommodating irregular cluster shapes and enabling parallel generation of posterior draws. To interpret the posterior over high-dimensional, noise-aware clusterings, we introduce a Wasserstein-type summarization method. Using a novel noise-aware loss, this approach compresses the large collection of posterior draws into a small set of weighted representatives, capturing distinct cluster structures and noise assignments. Further, we marginalise the uncertainty in the density level, avoiding reliance on a specific choice. We establish theoretical properties of the proposed loss and guarantees for the level-marginal posterior approximation. Results demonstrate the framework's effectiveness in capturing clustering and noise uncertainty.

Thyroid nodule malignancy assessment relies on ultrasound imaging, yet current risk stratification systems suffer from substantial interobserver variability, contributing to a high rate of unnecessary biopsies. We propose a radiomics-based framework combining quantitative feature extraction via PyRadiomics with a Generalized Additive Model (GAM) classifier, selected among six candidates for its favorable trade-off between discriminative performance and interpretability. To complement point predictions with uncertainty estimates, Cross-Conformal Prediction with label-conditional guarantees is applied, flagging ambiguous cases for clinical review rather than forcing a binary decision. Experiments are conducted on TN5000, a publicly available ultrasound dataset of thyroid nodules. After Spearman correlation filtering, 36 features are retained and the GAM achieves an AUC of $0.89$, a sensitivity of $0.83$, and a specificity of $0.82$ on the validation set. At a significance level $\epsilon=0.09$, selected to ensure empirical coverage guarantees on both classes, 63% of nodules on the test set receive a certain (i.e., single-label) prediction, with conditional precision of $0.94$ for certain malignant and conditional Negative Predictive Value (NPV) of $0.73$ for certain benign cases, while the remaining 37% are referred to the clinician.

Kernel Ridge Regression (KRR) is a theoretically sound algorithm for nonparametric learning, which is also often used to investigate properties of deep learning. However, the exact KRR formulation incurs time and memory costs preventing applicability to large datasets, motivating the development of approximate solvers. These algorithms rely on (i) iterative methods, to avoid the computational burden of direct solutions; (ii) sketching or batching, to reduce the problem size; and (iii) preconditioning, to accelerate convergence. While crucial for fast convergence, preconditioning incurs memory costs that can be managed by lowering the number of free parameters of the model, potentially at the expense of expressivity and hence accuracy. To address this challenge, we propose BLOKS, a Conjugate Gradient-based solver for Nyström-sketched KRR built around a novel preconditioner that combines a block Jacobi approximation with sketching. The BLOKS preconditioner has memory cost independent of both the dataset size and the number of parameters, allowing scaling to large models. We provide a memory-efficient, GPU-accelerated implementation and show that BLOKS achieves excellent accuracy while retaining favorable time/memory complexity. Indeed, on standard large-scale datasets used to benchmark kernel solvers, BLOKS achieves comparable or better performances than state-of-the-art methods.