Amortized Bandwidth Learning for Kernel Density Estimation under Logarithmic Score

# 1Introduction
Source: [https://arxiv.org/html/2608.20445](https://arxiv.org/html/2608.20445)
Amortized Bandwidth Learning for Kernel Density Estimation under Logarithmic Score

Junyi Liang1Hailiang Du2,1,3,\*

1School of Mathematics, East China University of Science and Technology, 200237, Shanghai, China

2Department of Mathematical Sciences; Institute of Hazard, Risk and Resilience, Durham University, DH1 3LE, Durham, United Kingdom

3Data Science Institute; Global School of Sustainability, The London School of Economics and Political Science, WC2A 2AE, London, United Kingdom

\*Corresponding author:[hailiang\.du@durham\.ac\.uk](mailto:[email protected])

Abstract

Kernel density estimation converts finite samples into probability densities, but its performance depends critically on bandwidth selection\. Classical selectors prescribe the sample\-to\-bandwidth rule analytically or asymptotically, or solve a new optimization for each sample\. An amortized framework is proposed that instead learns this mapping across a distribution of density\-estimation tasks by optimizing the logarithmic score\. A truncated\-and\-renormalized bounded\-support formulation enables stable learning across heterogeneous tasks, while affine standardization allows a selector trained on a single reference interval to transfer across bounded intervals\. Experiments under Gaussian sampling, a multi\-family benchmark, and randomized Gaussian\-mixture training show that the amortized selector consistently and substantially outperforms Silverman’s rule, the Sheather–Jones selector, and least\-squares cross\-validation, with especially large gains in small and heterogeneous samples\. Finite Gaussian mixtures provide a generic training mechanism supported by theirL1L^\{1\}approximation property\. Selectors trained in this way generalize strongly across different density structures, allowing the same trained selector to be applied directly to finite samples from unknown densities without specifying or fitting a distributional family\. This combination of broad applicability and strong empirical performance makes the framework attractive for a wide range of applications in which finite samples or ensembles must be converted into continuous probability densities\.

Keywords:amortized learning; bandwidth selection; Gaussian mixture models; kernel density estimation; logarithmic score

## 1Introduction

In many applications, one has access only to a finite sample or a finite ensemble from an underlying distribution, while the quantity ultimately needed for inference, evaluation, or decision\-making is a continuous probability density\. Examples include probabilistic forecasting, simulation\-based uncertainty quantification, and other settings in which tail probabilities, predictive likelihoods, or calibrated distributional summaries are of primary interest\. Kernel density estimation \(KDE\) is a classical nonparametric method for converting such finite samples into a density estimate\([18](https://arxiv.org/html/2608.20445#bib.bib1);[16](https://arxiv.org/html/2608.20445#bib.bib2)\)\. In kernel density estimation, the central practical difficulty is bandwidth selection, since the bandwidth controls the degree of smoothing and thereby has a decisive effect on estimation quality\. The present study focuses on the widely used Gaussian\-kernel setting\. Although the methodology is developed and evaluated in this setting, the underlying idea of learning a reusable bandwidth rule is not restricted to the Gaussian kernel and can be extended to other kernel families\.

A large literature has developed classical bandwidth selectors for kernel density estimation, including rule\-of\-thumb methods\([18](https://arxiv.org/html/2608.20445#bib.bib1);[16](https://arxiv.org/html/2608.20445#bib.bib2)\), plug\-in methods\([17](https://arxiv.org/html/2608.20445#bib.bib6);[19](https://arxiv.org/html/2608.20445#bib.bib7)\), and cross\-validation methods\([15](https://arxiv.org/html/2608.20445#bib.bib3);[9](https://arxiv.org/html/2608.20445#bib.bib4);[2](https://arxiv.org/html/2608.20445#bib.bib5)\)\. Rule\-of\-thumb selectors, such as Silverman’s rule, provide simple closed\-form prescriptions based on normal\-reference arguments together with global scale summaries of the sample\([18](https://arxiv.org/html/2608.20445#bib.bib1)\), but this simplicity comes at the cost of relying on a fixed reference structure that may be poorly matched to skewed, heavy\-tailed, contaminated, or multimodal densities\. Plug\-in methods are more flexible in that they typically begin from an asymptotic approximation to an estimation criterion and replace the resulting unknown population quantities, such as curvature\-related functionals, by data\-based estimates\([17](https://arxiv.org/html/2608.20445#bib.bib6);[19](https://arxiv.org/html/2608.20445#bib.bib7)\), but they still depend on asymptotic approximations and on the reliable estimation of higher\-order quantities, which may be unstable in small to moderate samples\. Cross\-validation methods rely less on an explicit reference model and instead select the bandwidth by optimizing a sample\-based criterion over candidate values\([15](https://arxiv.org/html/2608.20445#bib.bib3);[9](https://arxiv.org/html/2608.20445#bib.bib4);[2](https://arxiv.org/html/2608.20445#bib.bib5)\)\. These approaches differ in how the bandwidth is computed\. Rule\-of\-thumb and plug\-in procedures define reusable sample\-to\-bandwidth algorithms whose structure is prescribed analytically or asymptotically, whereas cross\-validation procedures typically solve a new optimization problem for each observed sample\. What they do not generally do is estimate the bandwidth\-selection mapping itself from a distribution of related density\-estimation tasks\. Consequently, they cannot directly exploit systematic cross\-task variation to learn which sample characteristics are predictive of good bandwidth choices under a specified evaluation criterion\. This motivates a different formulation in which the bandwidth rule is learned across tasks\. The proposed approach does not claim that classical selectors are unusable on new samples; rather, it replaces an analytically prescribed or repeatedly optimized selector with a mapping whose parameters are estimated from a collection of density\-estimation problems\. The resulting mapping remains adaptive to the observed sample, while also incorporating regularities identified across the training task distribution\.

Classical bandwidth selection in kernel density estimation has largely been developed under integrated squared\-error criteria, especially MISE and its asymptotic approximations\([18](https://arxiv.org/html/2608.20445#bib.bib1);[16](https://arxiv.org/html/2608.20445#bib.bib2);[19](https://arxiv.org/html/2608.20445#bib.bib7)\)\. These criteria are mathematically convenient, but KDE ultimately produces a probability density, so bandwidth quality is naturally viewed as a problem of probabilistic assessment\. Strict propriety alone does not determine how imperfect density estimates should be compared, since different strictly proper scoring rules can rank them differently\. Locality provides a further principle: for continuous densities, the logarithmic score is the only proper local scoring rule, depending directly on the density assigned to the realized outcome\([1](https://arxiv.org/html/2608.20445#bib.bib8);[3](https://arxiv.org/html/2608.20445#bib.bib9)\)\. By contrast, nonlocal scores can reward features of the reported density away from the realized outcome and may produce unfortunate evaluations\. The logarithmic score also has a direct interpretation in terms of probabilities and bits of information, while relative logarithmic score comparisons are invariant under smooth one\-to\-one transformations of the variable\([4](https://arxiv.org/html/2608.20445#bib.bib10)\)\. For these reasons, the proposed framework learns the bandwidth rule directly under the logarithmic score\.

These considerations suggest a different formulation of the KDE bandwidth problem\. Rather than fixing the functional form of the selector through an analytical or asymptotic prescription, or solving a new bandwidth\-optimization problem for every observed sample, the present work learns the sample\-to\-bandwidth mapping from a distribution of density\-estimation tasks\. The core idea is to learn a shared sample\-to\-bandwidth mapping that extracts structural information from each sample and outputs a bandwidth adapted to that task\. The resulting mapping, referred to as the amortized bandwidth selector, is reusable in the same operational sense as a classical selector, but its functional behaviour is estimated from cross\-task data rather than prescribed in advance\. Since the estimated object is a probability density, the amortized selector is trained and evaluated under the logarithmic score rather than an integrated squared\-error criterion\. To prevent the bandwidth from being driven to excessively large values by rare realizations in remote low\-density regions under logarithmic score evaluation, the framework adopts a bounded\-support formulation based on truncation and renormalization\. This restricts both the target density and the KDE to a common interval, so that the learning objective is concentrated on the region of practical interest rather than being influenced by behaviour outside it\. In this way, bandwidth selection is no longer treated as a fixed analytic prescription, but as a reusable rule learned directly from task variation\.

The resulting framework is evaluated from three complementary perspectives\. A Gaussian benchmark is first used to test whether the amortized selector captures the basic finite\-sample and scaling behaviour of KDE bandwidths in the most transparent setting\. A bounded multi\-family benchmark is then introduced to assess whether a single amortized selector remains effective under substantial structural variation across density families\. Finally, a Gaussian mixture model \(GMM\) benchmark examines whether a flexible Gaussian\-mixture generator can replace the designed multi\-family distribution as a source of tasks for learning the amortized selector\. Together, these experiments test whether the amortized selector improves empirical performance relative to classical selectors and remains effective across diverse density\-estimation tasks\.

The remainder of the paper is organized as follows\. Section 2 introduces the proposed bandwidth\-learning framework, including the amortized sample\-to\-bandwidth mapping, the logarithmic score objective, the bounded\-support formulation, and the interval\-standardization rule for cross\-interval transfer\. Section 3 presents the experimental results\. Section 4 concludes with a summary of the main findings and a discussion of limitations and future directions\.

## 2Bandwidth\-Learning Framework

The proposed framework retains KDE as the final density estimator and learns only the bandwidth\-selection rule\. Each observed sample is represented by a permutation\-invariant set of features\. A neural network maps these features to a bandwidth, which is then used to construct the KDE\. Unlike classical bandwidth selectors, the sample\-to\-bandwidth mapping is learned jointly across density\-estimation tasks and then reused for new samples without task\-specific reoptimization\.

The framework is developed in three parts\. Bandwidth selection is first formulated as an amortized learning problem, in which a shared predictor maps sample features to bandwidths across related tasks\. The training criterion is then specified as the logarithmic score, aligning the selector with the probabilistic objective used to assess density estimation\. Finally, a truncated\-and\-renormalized bounded\-support formulation is introduced together with affine standardization, which maps bounded task domains to a common reference interval and enables cross\-interval transfer through bandwidth rescaling\.

### 2\.1Amortized Bandwidth Prediction

This subsection formalizes the object that is learned\. Letffdenote an unknown density onℝ\\mathbb\{R\}, and let

S=\(x1,…,xn\)∼fnS=\(x\_\{1\},\\dots,x\_\{n\}\)\\sim f^\{n\}\(1\)be an observed i\.i\.d\. sample of sizenn\. For a bandwidthh\>0h\>0, the associated Gaussian\-kernel density estimator is

qh\(x;S\)=1n​h∑i=1nϕ\(x−xih\),ϕ\(u\)=\(2π\)−1/2e−u2/2\.q\_\{h\}\(x;S\)=\\frac\{1\}\{nh\}\\sum\_\{i=1\}^\{n\}\\phi\\\!\\left\(\\frac\{x\-x\_\{i\}\}\{h\}\\right\),\\qquad\\phi\(u\)=\(2\\pi\)^\{\-1/2\}e^\{\-u^\{2\}/2\}\.\(2\)The amortized selector treats the sample\-to\-bandwidth relationship as an unknown mapping to be learned across density\-estimation tasks\. Once trained, it can be applied to new samples without task\-specific reoptimization\. The bandwidth rule is written as

hθ​\(S\)=Fθ​\(feat⁡\(S\)\),h\_\{\\theta\}\(S\)=F\_\{\\theta\}\(\\mathrm\{feat\}\(S\)\),\(3\)wherefeat⁡\(S\)\\mathrm\{feat\}\(S\)is a finite\-dimensional vector of permutation\-invariant sample features,FθF\_\{\\theta\}is a parametric map with positive output, andθ\\thetadenotes its trainable parameters\. In the present proof\-of\-concept implementation,FθF\_\{\\theta\}is realized by a simple multilayer perceptron \(MLP\), and the feature map consists of the sample size, sample mean, sample standard deviation, sample skewness, and sample kurtosis:

feat⁡\(S\)=\(n,mean⁡\(S\),sd⁡\(S\),skew⁡\(S\),kurt⁡\(S\)\)\.\\mathrm\{feat\}\(S\)=\\bigl\(n,\\mathrm\{mean\}\(S\),\\mathrm\{sd\}\(S\),\\mathrm\{skew\}\(S\),\\mathrm\{kurt\}\(S\)\\bigr\)\.\(4\)These features capture basic information on sample size, location, scale, asymmetry, and tail behaviour while remaining low\-dimensional and interpretable\. For the Gaussian diagnostic benchmark below, a simplified special case is used in which the network learns only a dimensionless bandwidth ratio as a function of sample size\.

### 2\.2Logarithmic Score Criterion

The training criterion is taken from the framework of proper scoring rules\. A scoring rule assigns a numerical loss to a reported probability density and an observed outcome, with lower values indicating better probabilistic performance\. Strict propriety ensures that reporting the data\-generating distribution is optimal in expectation\([5](https://arxiv.org/html/2608.20445#bib.bib11);[6](https://arxiv.org/html/2608.20445#bib.bib12);[3](https://arxiv.org/html/2608.20445#bib.bib9)\)\. This property is essential, but it does not uniquely determine how imperfect density estimates should be compared: different strictly proper scoring rules can rank them differently\([4](https://arxiv.org/html/2608.20445#bib.bib10)\)\.

A further distinction is provided by locality\. A local scoring rule depends only on the density assigned to the realized outcome\. For continuous densities, the logarithmic score is, up to affine equivalence, the only proper local scoring rule\([1](https://arxiv.org/html/2608.20445#bib.bib8);[3](https://arxiv.org/html/2608.20445#bib.bib9);[10](https://arxiv.org/html/2608.20445#bib.bib13)\)\. This property is particularly relevant to bandwidth selection, since nonlocal scores can reward features of the estimated density away from the realized outcome and may therefore produce “unfortunate” evaluations\([4](https://arxiv.org/html/2608.20445#bib.bib10)\)\. They may also change their preferences under smooth transformations of the variable\. By contrast, the logarithmic score directly evaluates the density assigned to the realized value and has a direct interpretation in terms of probabilities and bits of information\([7](https://arxiv.org/html/2608.20445#bib.bib14);[14](https://arxiv.org/html/2608.20445#bib.bib15)\)\. Its sample average coincides with the negative log\-likelihood per observation, while its expectation is the cross\-entropy\.

These considerations motivate the use of the logarithmic score as the training and evaluation criterion for bandwidth selection\. For an independent realizationX∼fX\\sim f, the task\-specific logarithmic score criterion is

ℒ\(h;S,f\)=𝔼X∼f\[−log2qh\(X;S\)\]=−∫log2qh\(x;S\)f\(x\)dx\.\\mathcal\{L\}\(h;S,f\)=\\mathbb\{E\}\_\{X\\sim f\}\\left\[\-\\log\_\{2\}q\_\{h\}\(X;S\)\\right\]=\-\\int\\log\_\{2\}q\_\{h\}\(x;S\)f\(x\)\\,dx\.\(5\)To formalize learning across density\-estimation tasks, letΠ\\Pidenote the task distribution over\(f,n\)\(f,n\)\. The expected objective of the amortized selector over the task distribution is

𝒥Π​\(θ\)=𝔼\(f,n\)∼Π​𝔼S∼fn​\[ℒ⁡\(hθ​\(S\),S,f\)\]\.\\mathcal\{J\}\_\{\\Pi\}\(\\theta\)=\\mathbb\{E\}\_\{\(f,n\)\\sim\\Pi\}\\mathbb\{E\}\_\{S\\sim f^\{n\}\}\\left\[\\mathcal\{L\}\\bigl\(h\_\{\\theta\}\(S\);S,f\\bigr\)\\right\]\.\(6\)In practice, the expectation in𝒥Π​\(θ\)\\mathcal\{J\}\_\{\\Pi\}\(\\theta\)is approximated using a finite sample of training tasks:

𝒥^N​\(θ\)=1N​∑j=1Nℒ⁡\(hθ​\(Sj\),Sj,fj\),\(fj,nj\)∼Π,Sj∼fjnj\.\\widehat\{\\mathcal\{J\}\}\_\{N\}\(\\theta\)=\\frac\{1\}\{N\}\\sum\_\{j=1\}^\{N\}\\mathcal\{L\}\\bigl\(h\_\{\\theta\}\(S\_\{j\}\);S\_\{j\},f\_\{j\}\\bigr\),\\qquad\(f\_\{j\},n\_\{j\}\)\\sim\\Pi,\\quad S\_\{j\}\\sim f\_\{j\}^\{\\,n\_\{j\}\}\.\(7\)

### 2\.3Bounded\-Support Formulation and Cross\-Interval Transfer

A common bounded interval serves several purposes in the present learning framework\. First, under logarithmic score evaluation, rare realizations in remote low\-density regions can exert a disproportionate influence on the training objective and drive the learned bandwidth towards excessive smoothing\. Restricting the criterion to a region of practical interest limits this effect and focuses learning on the part of the distribution that is relevant for evaluation\. Second, affine standardization maps bounded task domains to a common reference interval\. This places sample features on comparable scales across tasks, so the predictor need not learn separately how to accommodate arbitrary shifts and changes of scale\. The bounded formulation also aligns naturally with applications in which the variable is physically constrained or only relevant over a finite range\. For such variables, a Gaussian\-kernel KDE may assign probability mass outside the admissible region; truncation and renormalization instead define the estimated density on the chosen interval\.

For any densityffwith positive mass on\[A,B\]\[A,B\], letf\[A,B\]​\(x\)=f⁡\(x\)​𝟏\[A,B\]​\(x\)/∫ABf⁡\(u\)​𝑑uf\_\{\[A,B\]\}\(x\)=f\(x\)\\mathbf\{1\}\_\{\[A,B\]\}\(x\)/\\int\_\{A\}^\{B\}f\(u\)\\,dudenote its truncated\-and\-renormalized version on that interval\. The Gaussian KDE is truncated and renormalized analogously:

q\[A,B\],h​\(x,S\)=qh​\(x,S\)​𝟏\[A,B\]​\(x\)∫ABqh​\(u,S\)​𝑑u\.q\_\{\[A,B\],h\}\(x;S\)=\\frac\{q\_\{h\}\(x;S\)\\mathbf\{1\}\_\{\[A,B\]\}\(x\)\}\{\\int\_\{A\}^\{B\}q\_\{h\}\(u;S\)\\,du\}\.\(8\)Thus both the target density and the KDE are rewritten as probability densities on the same bounded interval\. The corresponding bounded logarithmic score criterion is

ℒ\[A,B\]​\(h,S,f\)=𝔼X∼f\[A,B\]​\[−log2⁡q\[A,B\],h​\(X,S\)\]\.\\mathcal\{L\}\_\{\[A,B\]\}\(h;S,f\)=\\mathbb\{E\}\_\{X\\sim f\_\{\[A,B\]\}\}\\left\[\-\\log\_\{2\}q\_\{\[A,B\],h\}\(X;S\)\\right\]\.\(9\)In applications, the interval\[A,B\]\[A,B\]is treated as a working domain rather than as an estimate of the full mathematical support\. When physical or operational bounds are available, they provide a natural choice\. Otherwise,\[A,B\]\[A,B\]may be specified from representative historical or reference data so as to cover the range over which the density estimate is intended to be used, with suitable margins to avoid truncating non\-negligible probability mass\. A data\-driven alternative is to expand the observed sample range using the distribution\-free rank identity that, for an i\.i\.d\. continuous sample of sizenn, a new observation falls outside the current range with probability2/\(n\+1\)2/\(n\+1\)\. A reference model, such as the uniform distribution, can then be used to translate this probability into an interval expansion\.

To transfer the predictor between bounded intervals, a sample on\[A,B\]\[A,B\]is first mapped to the reference interval\[−1,1\]\[\-1,1\]by

T⁡\(x\)=x−ma,m=A\+B2,a=B−A2\.T\(x\)=\\frac\{x\-m\}\{a\},\\qquad m=\\frac\{A\+B\}\{2\},\\qquad a=\\frac\{B\-A\}\{2\}\.\(10\)
The predictor trained on\[−1,1\]\[\-1,1\]is then applied toT⁡\(S\)T\(S\), and the resulting bandwidth is rescaled to the original interval:

hθ,\[A,B\]​\(S\)=a​Fθ​\(feat⁡\(T⁡\(S\)\)\)\.h\_\{\\theta,\[A,B\]\}\(S\)=aF\_\{\\theta\}\\\!\\left\(\\mathrm\{feat\}\(T\(S\)\)\\right\)\.\(11\)Proposition 1 shows that this bandwidth rescaling is exactly compatible with the bounded logarithmic score criterion\.

###### Proposition 1\.

LetS∈\[A,B\]nS\\in\[A,B\]^\{n\}, letT⁡\(S\)∈\[−1,1\]nT\(S\)\\in\[\-1,1\]^\{n\}be its standardized sample, and letf~\\widetilde\{f\}denote the density ofY=\(X−m\)/aY=\(X\-m\)/aforX∼f\[A,B\]X\\sim f\_\{\[A,B\]\}\. Ifh\[A,B\]=a​h\[−1,1\]h\_\{\[A,B\]\}=a\\,h\_\{\[\-1,1\]\}, then

ℒ\[A,B\]​\(h\[A,B\],S,f\)=ℒ\[−1,1\]​\(h\[−1,1\],T⁡\(S\),f~\)\+log2⁡a\.\\mathcal\{L\}\_\{\[A,B\]\}\\bigl\(h\_\{\[A,B\]\};S,f\\bigr\)=\\mathcal\{L\}\_\{\[\-1,1\]\}\\bigl\(h\_\{\[\-1,1\]\};T\(S\),\\widetilde\{f\}\\bigr\)\+\\log\_\{2\}a\.\(12\)Consequently, the task\-specific minimizing bandwidths on the two intervals correspond under the scalingh\[A,B\]=a​h\[−1,1\]h\_\{\[A,B\]\}=a\\,h\_\{\[\-1,1\]\}\.

###### Proof\.

Under the change of variablesx=m\+a​yx=m\+ay, the target density, the KDE, and its truncation normalizer transform as

f\[A,B\]​\(x\)\\displaystyle f\_\{\[A,B\]\}\(x\)=a−1​f~​\(y\),\\displaystyle=a^\{\-1\}\\widetilde\{f\}\(y\),q\[A,B\],h\[A,B\]​\(x,S\)\\displaystyle q\_\{\[A,B\],h\_\{\[A,B\]\}\}\(x;S\)=a−1​q\[−1,1\],h\[−1,1\]​\(y,T⁡\(S\)\),\\displaystyle=a^\{\-1\}q\_\{\[\-1,1\],h\_\{\[\-1,1\]\}\}\(y;T\(S\)\),Z\[A,B\],h\[A,B\]​\(S\)\\displaystyle Z\_\{\[A,B\],h\_\{\[A,B\]\}\}\(S\)=Z\[−1,1\],h\[−1,1\]​\(T​\(S\)\),\\displaystyle=Z\_\{\[\-1,1\],h\_\{\[\-1,1\]\}\}\(T\(S\)\),whereZ\[A,B\],h​\(S\)=∫ABqh​\(u,S\)​𝑑uZ\_\{\[A,B\],h\}\(S\)=\\int\_\{A\}^\{B\}q\_\{h\}\(u;S\)\\,dudenotes the truncation normalizing constant\. Substituting these intoℒ\[A,B\]\\mathcal\{L\}\_\{\[A,B\]\}, the integration Jacobian cancels the scale factor off\[A,B\]f\_\{\[A,B\]\}, while the factora−1a^\{\-1\}inside the logarithm contributes the additive constantlog2⁡a\\log\_\{2\}a\. Sincelog2⁡a\\log\_\{2\}ais independent of the bandwidth, the minimizing bandwidths rescale by the factoraa\. ∎

## 3Experimental Demonstration

This section evaluates the proposed bandwidth\-learning framework from three complementary perspectives\. The first benchmark considers Gaussian sampling and examines finite\-sample bandwidth behaviour in the simplest setting\. The second considers a bounded multi\-family task distribution and tests whether a single shared predictor remains effective across different distributional shapes\. The third uses randomized finite Gaussian mixtures as the training task distribution and examines whether a single flexible density generator can provide a broad source of tasks for learning a reusable bandwidth rule\.

Implementation choices are kept simple throughout\. The aim is to evaluate the bandwidth\-learning formulation itself rather than to optimize network architecture or numerical implementation\. Accordingly, the predictors use low\-dimensional sample features, with positivity of the bandwidth enforced through a softplus output transformation\. All logarithms use base 2, so logarithmic score differences have the usual interpretation in bits of information\. A one\-bit lower score corresponds to assigning, on average, twice as much probability density to the realized outcomes\([4](https://arxiv.org/html/2608.20445#bib.bib10)\)\.

### 3\.1Comparison of Bandwidth Selectors under Gaussian Sampling

This experiment compares the amortized bandwidth selector with three classical selectors under Gaussian sampling\. This controlled setting removes variation in skewness, tail behaviour, and modality, allowing the comparison to isolate the effect of sample size on logarithmic score performance\. Performance under non\-Gaussian density shapes is examined in the subsequent experiments\.

For this Gaussian diagnostic benchmark, the amortized selector is restricted to the scale\-equivariant form

hθ​\(S\)=sd⁡\(S\)​Fθ​\(n\),Fθ​\(n\)\>0\.h\_\{\\theta\}\(S\)=\\operatorname\{sd\}\(S\)F\_\{\\theta\}\(n\),\\qquad F\_\{\\theta\}\(n\)\>0\.\(13\)Thus, the network learns only a dimensionless bandwidth ratio as a function of sample size, while the underlying distributional shape is fixed\. The amortized selector is compared with Silverman’s rule, the Sheather–Jones selector, and least\-squares cross\-validation\. Their definitions and implementation details are provided in Appendix[A](https://arxiv.org/html/2608.20445#A1)\. All four selectors are translation invariant and scale equivariant, up to negligible numerical implementation error\. The corresponding derivations are given in Appendix[B](https://arxiv.org/html/2608.20445#A2)\. It is therefore sufficient to conduct the comparison under standard\-normal sampling, since translating or rescaling a Gaussian distribution does not alter the ordering of the selectors or their pairwise logarithmic score differences\.

Training samples are generated from𝒩⁡\(0,1\)\\mathcal\{N\}\(0,1\), withn∼Unif⁡\{5,…,256\}n\\sim\\operatorname\{Unif\}\\\{5,\\ldots,256\\\}\. Both the target density and the KDE are defined on the full real line; no truncation or renormalization is applied in this experiment\. The predictor is trained on online\-generated Gaussian tasks, with model selection based on a fixed validation set\. Details of the network architecture and training procedure are provided in Appendix[C](https://arxiv.org/html/2608.20445#A3)\. For each displayed sample size, performance is evaluated onNtest=30000N\_\{\\mathrm\{test\}\}=30000independent test samples\.

![Refer to caption](https://arxiv.org/html/2608.20445v1/fig_gaussian_amortized_selector_comparison.png)Figure 1:Empirical logarithmic score versus sample size under Gaussian sampling\. Scores are averaged overNtest=30000N\_\{\\mathrm\{test\}\}=30000independent samples at eachnn; shaded bands show 90% bootstrap intervals, and the black dashed line denotes the score of the true𝒩⁡\(0,1\)\\mathcal\{N\}\(0,1\)density\. The horizontal axis is linear forn≤20n\\leq 20and logarithmically compressed thereafter, while the vertical axis is linear up to4\.54\.5bits and logarithmically compressed above\.Figure[1](https://arxiv.org/html/2608.20445#S3.F1)shows that the amortized selector attains the lowest empirical logarithmic score throughout the displayed sample\-size range and remains consistently closer to the standard\-normal reference than the three classical selectors\. The poor performance of the classical selectors in very small samples is consistent with occasional severe undersmoothing\. Accidental clustering of observations can lead to an extremely small selected bandwidth, causing the resulting KDE to assign very low density to outcomes that remain plausible under the standard\-normal distribution and hence to incur large logarithmic penalties\. The amortized selector is less sensitive to accidental local spacing irregularities because its bandwidth ratio is learned across training tasks rather than determined by sample\-specific spacing patterns\. Training under the logarithmic score also discourages bandwidths that assign negligible density to outcomes that remain plausible under the data\-generating distribution\. This strong penalization is not an artefact of the logarithmic score, but an appropriate response to severe density\-estimation errors\([4](https://arxiv.org/html/2608.20445#bib.bib10)\)\. As the sample size increases, all four selectors move towards the true standard\-normal reference and their differences gradually narrow, although the amortized selector retains a clear advantage through at leastn=50n=50\.

### 3\.2Comparison of Bandwidth Selectors under Bounded Multi\-Family Sampling

The Gaussian comparison isolates finite\-sample bandwidth behaviour under a single regular density family\. The present experiment evaluates whether a shared amortized selector remains effective across a much broader collection of density shapes, including skewed, heavy\-tailed, boundary\-concentrated, multimodal, and multi\-scale distributions\.

The benchmark comprises ten equally weighted distribution families: Gaussian, Laplace, Student\-tt, Gamma, Beta, Logistic, Lognormal, Bimodal, Trimodal, and Spike–slab\. Family\-specific parameter ranges are chosen so that at least90%90\\%of the untruncated probability mass lies in\[−1,1\]\[\-1,1\], providing a practical balance between concentrating most probability mass within the working interval and retaining some tail mass outside it\. Similar results were obtained with alternative thresholds\. Further details are given in Table[1](https://arxiv.org/html/2608.20445#S3.T1)\. The bounded\-support formulation of Section 2\.3 is then applied\.

Table 1:Densities and parameter specifications for the bounded multi\-family benchmark\.FamilyDensityg⁡\(x\)g\(x\)ParametersGaussiang⁡\(x\)=ϕσ​\(x−μ\)g\(x\)=\\phi\_\{\\sigma\}\(x\-\\mu\)σ∈\[0\.13,0\.61\]\\sigma\\in\[0\.13,0\.61\]\.Laplaceg\(x\)=\(2b\)−1exp\{−\|x−μ\|/b\}g\(x\)=\(2b\)^\{\-1\}\\exp\\\{\-\|x\-\\mu\|/b\\\}b∈\[0\.06,0\.43\]b\\in\[0\.06,0\.43\]\.Student\-ttg⁡\(x\)=s−1​tν​\{\(x−μ\)/s\}g\(x\)=s^\{\-1\}t\_\{\\nu\}\\\{\(x\-\\mu\)/s\\\}ν∼L​U​\[3,30\]\\nu\\sim LU\[3,30\];s∈\[0\.02,0\.59\]s\\in\[0\.02,0\.59\]\.Gammay=k​θ\+ε⁡\(x−μ\),g\(x\)=γk,θ\(y\)𝟏\{y\>0\}\.\\begin\{array\}\[\]\{@\{\}l@\{\}\}y=k\\theta\+\\varepsilon\(x\-\\mu\),\\\\\[\-2\.0pt\] g\(x\)=\\gamma\_\{k,\\theta\}\(y\)\\mathbf\{1\}\_\{\\\{y\>0\\\}\}\.\\end\{array\}k∼L​U​\[0\.5,10\]k\\sim LU\[0\.5,10\];ε∈\{−1,1\}\\varepsilon\\in\\\{\-1,1\\\}equiprobably;θ∈\[0\.04,1\.4\]\\theta\\in\[0\.04,1\.4\]\.Betaz=aa\+b\+x−μs,g\(x\)=s−1βa,b\(z\)𝟏\{0<z<1\}\.\\begin\{array\}\[\]\{@\{\}l@\{\}\}z=\\dfrac\{a\}\{a\+b\}\+\\dfrac\{x\-\\mu\}\{s\},\\\\\[\-2\.0pt\] g\(x\)=s^\{\-1\}\\beta\_\{a,b\}\(z\)\\mathbf\{1\}\_\{\\\{0<z<1\\\}\}\.\\end\{array\}Four equiprobable regimes:a,b∼iidU⁡\[0\.4,0\.9\]a,b\\stackrel\{\{\\scriptstyle\\mathrm\{iid\}\}\}\{\{\\sim\}\}U\[0\.4,0\.9\],U⁡\[1\.2,5\]U\[1\.2,5\], orU⁡\[0\.8,1\.2\]U\[0\.8,1\.2\]; oraaandbbare drawn, in random order, fromU⁡\[0\.4,0\.9\]U\[0\.4,0\.9\]andU⁡\[1\.2,5\]U\[1\.2,5\]\.s∈\[0\.7,8\.3\]s\\in\[0\.7,8\.3\]\.Logisticg⁡\(x\)=exp\{−\(x−μ\)/s\}s\[1\+exp\{−\(x−μ\)/s\}\]2g\(x\)=\\frac\{\\exp\\\{\-\(x\-\\mu\)/s\\\}\}\{s\[1\+\\exp\\\{\-\(x\-\\mu\)/s\\\}\]^\{2\}\}s∈\[0\.06,0\.34\]s\\in\[0\.06,0\.34\]\.Lognormaly=eτ2/2\+ε⁡\(x−μ\)/s,g\(x\)=s−1ℓτ\(y\)𝟏\{y\>0\}\.\\begin\{array\}\[\]\{@\{\}l@\{\}\}y=e^\{\\tau^\{2\}/2\}\+\\varepsilon\(x\-\\mu\)/s,\\\\\[\-2\.0pt\] g\(x\)=s^\{\-1\}\\ell\_\{\\tau\}\(y\)\\mathbf\{1\}\_\{\\\{y\>0\\\}\}\.\\end\{array\}τ∼U⁡\[0\.25,1\]\\tau\\sim U\[0\.25,1\];ε∈\{−1,1\}\\varepsilon\\in\\\{\-1,1\\\}equiprobably;s∈\[0\.01,2\.4\]s\\in\[0\.01,2\.4\]\.Bimodalg⁡\(x\)=w​ϕσ1​\(x−μ1\)\+\(1−w\)​ϕσ2​\(x−μ2\)g\(x\)=w\\phi\_\{\\sigma\_\{1\}\}\(x\-\\mu\_\{1\}\)\+\(1\-w\)\\phi\_\{\\sigma\_\{2\}\}\(x\-\\mu\_\{2\}\)w∼U⁡\[0\.2,0\.8\]w\\sim U\[0\.2,0\.8\]\.\(m1,m2\)=\(−1,1\),m¯=w​m1\+\(1−w\)​m2\\begin\{array\}\[\]\{@\{\}l@\{\}\}\(m\_\{1\},m\_\{2\}\)=\(\-1,1\),\\\\\[\-1\.0pt\] \\bar\{m\}=wm\_\{1\}\+\(1\-w\)m\_\{2\}\\end\{array\}ρ1,ρ2∼iidL​U​\[0\.15,0\.45\]\\rho\_\{1\},\\rho\_\{2\}\\stackrel\{\{\\scriptstyle\\mathrm\{iid\}\}\}\{\{\\sim\}\}LU\[0\.15,0\.45\]\.μj=μ\+λ⁡\(mj−m¯\),σj=λ​ρj\\begin\{array\}\[\]\{@\{\}l@\{\}\}\\mu\_\{j\}=\\mu\+\\lambda\(m\_\{j\}\-\\bar\{m\}\),\\\\\[\-1\.0pt\] \\sigma\_\{j\}=\\lambda\\rho\_\{j\}\\end\{array\}μ1∈\[−1\.0,0\.4\]\\mu\_\{1\}\\in\[\-1\.0,0\.4\],μ2∈\[−0\.4,1\.0\]\\mu\_\{2\}\\in\[\-0\.4,1\.0\];σ1,σ2∈\[0\.04,0\.34\]\\sigma\_\{1\},\\sigma\_\{2\}\\in\[0\.04,0\.34\]\.Trimodalg⁡\(x\)=∑j=13wj​ϕσj​\(x−μj\)g\(x\)=\\sum\_\{j=1\}^\{3\}w\_\{j\}\\phi\_\{\\sigma\_\{j\}\}\(x\-\\mu\_\{j\}\)\(w1,w2,w3\)∼Dirichlet⁡\(2,2,2\)\(w\_\{1\},w\_\{2\},w\_\{3\}\)\\sim\\operatorname\{Dirichlet\}\(2,2,2\),minj⁡wj≥0\.1\\min\_\{j\}w\_\{j\}\\geq 0\.1\.\(m1,m2,m3\)=\(−dL,0,dR\)\(m\_\{1\},m\_\{2\},m\_\{3\}\)=\(\-d\_\{L\},0,d\_\{R\}\)dL,dR∼iidU⁡\[0\.8,1\.2\]d\_\{L\},d\_\{R\}\\stackrel\{\{\\scriptstyle\\mathrm\{iid\}\}\}\{\{\\sim\}\}U\[0\.8,1\.2\]\.m¯=∑j=13wj​mj\\bar\{m\}=\\sum\_\{j=1\}^\{3\}w\_\{j\}m\_\{j\}ρj∼iidL​U​\[0\.1,0\.3\]\\rho\_\{j\}\\stackrel\{\{\\scriptstyle\\mathrm\{iid\}\}\}\{\{\\sim\}\}LU\[0\.1,0\.3\]\.μj=μ\+λ⁡\(mj−m¯\),σj=λ​ρj\\begin\{array\}\[\]\{@\{\}l@\{\}\}\\mu\_\{j\}=\\mu\+\\lambda\(m\_\{j\}\-\\bar\{m\}\),\\\\\[\-1\.0pt\] \\sigma\_\{j\}=\\lambda\\rho\_\{j\}\\end\{array\}μ1∈\[−1\.3,0\.4\]\\mu\_\{1\}\\in\[\-1\.3,0\.4\],μ2∈\[−0\.7,0\.7\]\\mu\_\{2\}\\in\[\-0\.7,0\.7\],μ3∈\[−0\.4,1\.4\]\\mu\_\{3\}\\in\[\-0\.4,1\.4\];σ1∈\[0\.03,0\.29\]\\sigma\_\{1\}\\in\[0\.03,0\.29\],σ2∈\[0\.03,0\.31\]\\sigma\_\{2\}\\in\[0\.03,0\.31\],σ3∈\[0\.03,0\.28\]\\sigma\_\{3\}\\in\[0\.03,0\.28\]\.Spike–slabg⁡\(x\)=w​ϕσspike​\(x−μ\)\+\(1−w\)​ϕσslab​\(x−μ\)\\begin\{array\}\[\]\{@\{\}l@\{\}\}g\(x\)=w\\phi\_\{\\sigma\_\{\\mathrm\{spike\}\}\}\(x\-\\mu\)\\\\\[\-2\.0pt\] \\qquad\+\(1\-w\)\\phi\_\{\\sigma\_\{\\mathrm\{slab\}\}\}\(x\-\\mu\)\\end\{array\}w∼U⁡\[0\.2,0\.8\]w\\sim U\[0\.2,0\.8\]\.\(σspike,σslab\)=λ⁡\(ρ,1\)\(\\sigma\_\{\\mathrm\{spike\}\},\\sigma\_\{\\mathrm\{slab\}\}\)=\\lambda\(\\rho,1\)ρ∼L​U​\[0\.05,0\.20\]\\rho\\sim LU\[0\.05,0\.20\];σspike∈\[0\.01,0\.26\]\\sigma\_\{\\mathrm\{spike\}\}\\in\[0\.01,0\.26\];σslab∈\[0\.12,1\.43\]\\sigma\_\{\\mathrm\{slab\}\}\\in\[0\.12,1\.43\]\.Note:Each family contributes10,00010\{,\}000independently generated settings \(100,000100\{,\}000in total\)\. For each setting, independentμ∼U⁡\[−0\.5,0\.5\]\\mu\\sim U\[\-0\.5,0\.5\]andc∼U⁡\[0\.9,1\]c\\sim U\[0\.9,1\]are drawn\. Given the family\-specific shape parameters, the overall scale \(λ\\lambdafor mixtures\) is chosen numerically to satisfy∫−11g⁡\(u\)​𝑑u=c\\int\_\{\-1\}^\{1\}g\(u\)\\,du=c\. Intervals without an explicit sampling distribution are realised ranges over the10,00010\{,\}000settings for that family\. Sampling and evaluation usef\[−1,1\]​\(x\)=g⁡\(x\)​𝟏\[−1,1\]​\(x\)/∫−11g⁡\(u\)​𝑑uf\_\{\[\-1,1\]\}\(x\)=g\(x\)\\mathbf\{1\}\_\{\[\-1,1\]\}\(x\)/\\int\_\{\-1\}^\{1\}g\(u\)\\,du\. Here,ϕσ\\phi\_\{\\sigma\},tνt\_\{\\nu\},γk,θ\\gamma\_\{k,\\theta\},βa,b\\beta\_\{a,b\}, andℓτ\\ell\_\{\\tau\}denote the Gaussian, standard Student\-tt, Gamma, Beta, and lognormal densities;U⁡\[a,b\]≡Unif⁡\[a,b\]U\[a,b\]\\equiv\\operatorname\{Unif\}\[a,b\]andL​U​\[a,b\]≡LogUnif⁡\[a,b\]LU\[a,b\]\\equiv\\operatorname\{LogUnif\}\[a,b\]\. Displayed realised bounds are rounded to one or two decimal places; all computations use full\-precision values\.For each family, 10,000 parameter settings are generated in advance and sampled during training, withn∼Unif⁡\{5,…,256\}n\\sim\\operatorname\{Unif\}\\\{5,\\ldots,256\\\}\. The amortized selector uses the feature representation in \(4\); further implementation details are provided in Appendix[C](https://arxiv.org/html/2608.20445#A3)\. All comparisons are conducted on\[−1,1\]\[\-1,1\], since Proposition 1 together with the scale equivariance of the selectors implies that the results transfer directly to any affine image of this interval\. At each displayed sample size, performance is averaged equally across the ten families using 3000 independently generated test tasks per family\.

![Refer to caption](https://arxiv.org/html/2608.20445v1/fig_multifamily_amortized_selector_comparison_90coverage.png)Figure 2:Empirical logarithmic score versus sample size under bounded multi\-family sampling on\[−1,1\]\[\-1,1\]\. Scores are averaged equally across the ten families usingNtest=3000N\_\{\\mathrm\{test\}\}=3000independent tasks per family at each sample size\. Shaded bands show 90% bootstrap intervals, and the black line denotes the corresponding score of the true target densities\.Figure[2](https://arxiv.org/html/2608.20445#S3.F2)shows that the amortized selector attains the lowest empirical logarithmic score throughout the displayed sample\-size range and remains consistently closer to the true\-density reference than the three classical selectors\. Its advantage is more pronounced than under Gaussian sampling, indicating that the improvement extends to heterogeneous density shapes\. The stronger separation is consistent with the fact that the classical selectors are neither adapted to this task distribution nor optimized directly under the logarithmic score\.

### 3\.3Comparison of Bandwidth Selectors under Bounded GMM Sampling

The multi\-family benchmark shows that a single amortized selector can remain effective across a broad range of density families\. Its construction, however, requires an explicit choice of which distribution families to include\. This motivates examining whether a single flexible density generator can provide a similarly broad training distribution\.

Gaussian mixture models provide a natural candidate for this purpose\. A finite Gaussian mixture has the form

gK​\(x\)=∑k=1Kwk​ϕσk​\(x−μk\),wk≥0,∑k=1Kwk=1\.g\_\{K\}\(x\)=\\sum\_\{k=1\}^\{K\}w\_\{k\}\\phi\_\{\\sigma\_\{k\}\}\(x\-\\mu\_\{k\}\),\\qquad w\_\{k\}\\geq 0,\\qquad\\sum\_\{k=1\}^\{K\}w\_\{k\}=1\.\(14\)[12](https://arxiv.org/html/2608.20445#bib.bib16)show that finite location–scale mixtures generated from an essentially bounded base density can approximate target densities arbitrarily well inLpL^\{p\}\. Since the Gaussian density satisfies this condition, finite Gaussian mixtures can approximate arbitrary probability densities inL1L^\{1\}as the number and parameters of the components vary\. This provides a theoretical motivation for using GMMs as a flexible source of density\-estimation tasks\. This approximation property is preserved under truncation and renormalization on\[−1,1\]\[\-1,1\]\.

Randomized finite Gaussian mixtures are generated for three component counts,K∈\{8,16,32\}K\\in\\\{8,16,32\\\}\. For each fixedKK, the mixture weights and component locations are generated according to

α\\displaystyle\\alpha∼LogUnif⁡\[0\.5,4\],\\displaystyle\\sim\\operatorname\{LogUnif\}\[0\.5,4\],\(15\)\(w1,…,wK\)\\displaystyle\(w\_\{1\},\\ldots,w\_\{K\}\)∼Dirichlet⁡\(αK,…,αK\),\\displaystyle\\sim\\operatorname\{Dirichlet\}\\left\(\\frac\{\\alpha\}\{K\},\\ldots,\\frac\{\\alpha\}\{K\}\\right\),μk\\displaystyle\\mu\_\{k\}∼iidUnif⁡\[−1,1\]\.\\displaystyle\\stackrel\{\{\\scriptstyle\\mathrm\{iid\}\}\}\{\{\\sim\}\}\\operatorname\{Unif\}\[\-1,1\]\.The use ofα/K\\alpha/Kkeeps the total Dirichlet concentration equal toα\\alphafor allKK, so that changingKKdoes not simultaneously alter the overall concentration of the mixture weights\. This scaling also has a theoretical motivation: finite mixtures with symmetric Dirichlet weights of this form approach Dirichlet\-process mixtures with concentration parameterα\\alphaasKKincreases\([8](https://arxiv.org/html/2608.20445#bib.bib17)\)\. Thus,KKcontrols the number of available mixture components, whereasα\\alphacontrols the degree of weight imbalance across tasks\. The rangeα∈\[0\.5,4\]\\alpha\\in\[0\.5,4\], sampled on the logarithmic scale, is chosen to generate substantial variation in mixture\-weight concentration rather than to represent theoretically optimal values\. The component locations are sampled independently over the reference interval without imposing separation or regular\-spacing constraints, allowing components to overlap naturally so that the nominal component countKKdoes not directly determine the number of modes\.

Component scales are generated in two stages\. A task\-level reference scale and a within\-task log\-scale dispersion are first sampled as

σ0∼LogUnif⁡\[0\.025,0\.22\],τ∼Unif⁡\[0\.10,1\.00\]\.\\sigma\_\{0\}\\sim\\operatorname\{LogUnif\}\[0\.025,0\.22\],\\qquad\\tau\\sim\\operatorname\{Unif\}\[0\.10,1\.00\]\.\(16\)Conditional on\(σ0,τ\)\(\\sigma\_\{0\},\\tau\), the component log\-scales are sampled according to

log⁡σk∼iid𝒩⁡\(log⁡σ0,τ2\),\\log\\sigma\_\{k\}\\stackrel\{\{\\scriptstyle\\mathrm\{iid\}\}\}\{\{\\sim\}\}\\mathcal\{N\}\\left\(\\log\\sigma\_\{0\},\\tau^\{2\}\\right\),\(17\)with sampling restricted toσk∈\[0\.015,0\.45\]\\sigma\_\{k\}\\in\[0\.015,0\.45\]\. The reference scaleσ0\\sigma\_\{0\}determines the typical component width within a task, whereasτ\\taucontrols the heterogeneity of component widths\. Together, they allow both overall and local scales to vary across tasks\. The bounds onσk\\sigma\_\{k\}exclude extremely narrow or diffuse components relative to the reference interval\[−1,1\]\[\-1,1\]\. These ranges are chosen to provide broad but nondegenerate variation\. The bounded\-support formulation of Section 2\.3 is then applied to each generated mixture on\[−1,1\]\[\-1,1\]\.

The GMM construction used here is not tuned to the test distributions and is not intended to be optimal\. Empirically, similar performance was obtained under reasonable variations of the GMM construction\. Other common GMM constructions are designed primarily for mixture inference rather than for generating a heterogeneous collection of density\-estimation tasks\. Bayesian mixture models, for example, often introduce hierarchical or shrinkage priors to regularize component parameters or to encourage a smaller number of occupied components\([13](https://arxiv.org/html/2608.20445#bib.bib18);[11](https://arxiv.org/html/2608.20445#bib.bib19)\)\. Such constructions are well suited to mixture inference but impose additional structure on component weights or scales, while equal\-weight mixtures eliminate variation in the relative importance of the components\. The present construction instead uses direct randomization to allow weight imbalance, component overlap, overall scale, and within\-task scale heterogeneity to vary across tasks\.

Because overlapping components can produce fewer modes than the nominal component countKK, a Monte Carlo diagnostic is conducted using50005000independently generated densities for each value ofKK\. Table[2](https://arxiv.org/html/2608.20445#S3.T2)reports the empirical distribution of the detected number of modes\. The realized modal complexity changes only moderately withKK\. Across all three generators, two or three modes are most common, and most densities contain no more than four modes\. Even forK=32K=32, only6\.0%6\.0\\%of densities have six or more modes\. Thus, increasingKKadds available mixture components without producing a proportional increase in modal complexity\.

Table 2:Empirical distributions of detected mode counts under the bounded GMM generators\. Entries are percentages from50005000independently generated densities for eachKK\.Number of detected modesNumber of componentsKK12345678≥9\\geq 9817\.0%36\.2%27\.2%13\.7%4\.8%1\.0%0\.2%0\.0%0\.0%1614\.4%32\.0%26\.8%15\.2%7\.3%3\.1%0\.9%0\.2%0\.0%3213\.5%30\.0%26\.3%15\.5%8\.7%4\.0%1\.2%0\.6%0\.2%Note:Modes are detected numerically using a2%2\\%relative\-prominence threshold; boundary modes are included\.Three amortized selectors are then trained separately using theK=8K=8,K=16K=16, andK=32K=32GMM generators, with a fresh GMM generated independently for each training task\. Sample sizes are drawn uniformly from\{5,…,256\}\\\{5,\\ldots,256\\\}, and the selector uses the same feature representation, bounded logarithmic score objective, network architecture, and model\-selection procedure as in the multi\-family experiment\. Implementation details are provided in Appendix[C](https://arxiv.org/html/2608.20445#A3)\.

Evaluation follows the protocol of Section[3\.2](https://arxiv.org/html/2608.20445#S3.SS2)\. The three GMM\-trained selectors are tested on independently generated boundedK=32K=32GMM tasks and on the same bounded multi\-family benchmark\. The previously trained multi\-family amortized selector is included in both comparisons\.

![Refer to caption](https://arxiv.org/html/2608.20445v1/fig_gmmk32_amortized_selector_comparison.png)Figure 3:Empirical logarithmic score versus sample size under boundedK=32K=32GMM sampling on\[−1,1\]\[\-1,1\]\. Scores are averaged overNtest=30000N\_\{\\mathrm\{test\}\}=30000independent GMM tasks at each sample size\. GMMK=8K=8, GMMK=16K=16, and GMMK=32K=32denote amortized selectors trained using the respective GMM generators, while Multi\-family denotes the selector trained on the bounded multi\-family task distribution\. Shaded bands show90%90\\%bootstrap intervals, and the black line denotes the corresponding score of the true target densities\.Figure[3](https://arxiv.org/html/2608.20445#S3.F3)shows a more pronounced separation between the GMM\-trained amortized selectors and the classical selectors than in the bounded multi\-family experiment, particularly at small and moderate sample sizes\. All three GMM\-trained selectors remain substantially closer to the true\-density reference throughout the displayed sample\-size range, with the gap narrowing asnnincreases\. The multi\-family\-trained selector also performs well on the GMM test distribution, although its empirical logarithmic score lies slightly above the GMM\-trained selectors\. This provides further evidence of cross\-distribution generalization\.

![Refer to caption](https://arxiv.org/html/2608.20445v1/fig_multifamily_gmm_amortized_selector_comparison.png)Figure 4:Empirical logarithmic score versus sample size under bounded multi\-family sampling on\[−1,1\]\[\-1,1\]\. Scores are averaged equally across the ten families usingNtest=3000N\_\{\\mathrm\{test\}\}=3000independent tasks per family at each sample size\. GMMK=8K=8,1616, and3232denote selectors trained using the corresponding GMM generators, while Multi\-family denotes the multi\-family\-trained selector\. Shaded bands show90%90\\%bootstrap intervals, and the black line denotes the corresponding score of the true target densities\.Figure[4](https://arxiv.org/html/2608.20445#S3.F4)shows that the selector trained directly on the multi\-family task distribution attains the lowest empirical logarithmic score among the learned selectors on its matched test distribution\. Nevertheless, all three GMM\-trained selectors remain close to the multi\-family\-trained selector and substantially outperform the classical bandwidth selectors, particularly at small and moderate sample sizes\. These results show that the GMM\-trained selectors generalize effectively to the multi\-family benchmark, even though the ten distribution families are not explicitly included in their training task distributions\.

Across both the GMM and multi\-family test distributions, the three GMM\-trained selectors produce very similar empirical logarithmic scores, indicating that performance is not strongly sensitive to the component count over the range considered\. The mode\-count diagnostic in Table[2](https://arxiv.org/html/2608.20445#S3.T2)provides complementary evidence that the realized modal complexity changes only moderately asKKincreases\.

Together, theL1L^\{1\}approximation property and the empirical cross\-distribution results support randomized finite GMMs as a flexible training distribution for amortized bandwidth learning\.

## 4Summary and Discussion

This paper develops an amortized approach to bandwidth selection for Gaussian\-kernel density estimation\. Classical bandwidth selectors obtain the bandwidth through analytically or asymptotically prescribed rules or through repeated sample\-specific optimization, as represented here by Silverman’s rule, the Sheather–Jones plug\-in selector, and least\-squares cross\-validation\. In contrast, the proposed framework learns the sample\-to\-bandwidth mapping itself across a distribution of density\-estimation tasks\. The density estimator remains classical KDE; only the bandwidth\-selection rule is learned\. Training is carried out directly under the logarithmic score, aligning bandwidth selection with the probabilistic criterion used to evaluate the resulting density estimate\.

A key part of the framework is the bounded\-support formulation\. Without restricting training to a common region of interest, rare realizations in remote low\-density regions can dominate the logarithmic score and drive the learned bandwidth towards excessive smoothing, making it difficult to train a single selector effectively across heterogeneous tasks\. The truncated\-and\-renormalized formulation addresses this problem by placing the target density and the KDE on a common bounded interval and concentrating the learning objective on the region of practical interest\. This restriction is not particularly limiting in practice\. The bounded interval is treated as a working domain rather than as an estimate of the full mathematical support\. Physical or operational bounds provide a natural choice when available; otherwise, the interval may be specified from representative historical or reference data with suitable margins, or constructed from a data\-driven expansion of the observed sample range\. Proposition 1 further shows that affine standardization changes the bounded logarithmic score only by an additive constant and rescales the corresponding minimizing bandwidth, allowing a selector trained on the reference interval to be transferred by standardization and bandwidth rescaling\.

The empirical results show a consistent advantage of amortized bandwidth selection over the classical selectors considered here\. Under Gaussian sampling, a setting particularly favourable to classical bandwidth selection and especially to Silverman’s normal\-reference rule, the amortized selector nevertheless attains the lowest empirical logarithmic score throughout the examined sample\-size range, with the largest gains at small sample sizes\. The multi\-family experiment further strengthens this result\. A single amortized selector maintains strong performance across widely varying density structures and shows an even larger advantage over the classical selectors\.

The GMM experiment further shows that a manually enumerated collection of distribution families is not required to construct an effective training task distribution\. Finite Gaussian mixtures provide a flexible density\-approximation framework, supported by theirL1L^\{1\}approximation property, and allow a broad range of density structures to be generated within a single generative framework\. The GMM\-trained selectors show little sensitivity to the component count over the range considered and substantially outperform the classical selectors on independently generated GMM tasks\. More importantly, they remain close to the matched multi\-family\-trained selector on the same benchmark, even though those distribution families are not explicitly included in the GMM training distribution\. Conversely, the multi\-family\-trained selector also performs well on the GMM benchmark\. Once trained, the same amortized selector can be applied directly to finite samples from unknown densities without specifying, identifying, or fitting their distributional families\. The cross\-distribution experiments show that this broad applicability is accompanied by strong empirical performance across diverse density structures\.

The significance of amortization is therefore not merely computational\. Once trained, the selector requires only a single forward evaluation for each new sample, but the more fundamental change is statistical: information accumulated across a distribution of density\-estimation tasks is used to learn the bandwidth\-selection mapping itself\. This allows recurring relationships between observable sample characteristics and appropriate smoothing to be exploited rather than rediscovered from each finite sample in isolation\. The especially large gains in small and heterogeneous samples are consistent with this interpretation, since these are precisely the settings in which an individual sample contains limited information about the underlying density structure\. An interactive implementation of the amortized bandwidth selector is available through the[online application](https://amortized-kde.streamlit.app/), where the trained selector can be applied directly to user\-supplied one\-dimensional samples\.

The present study is restricted to one\-dimensional Gaussian\-kernel KDE with a single global bandwidth and a small collection of interpretable sample features\. These choices make it possible to isolate the effect of learning the bandwidth rule itself, but they also leave substantial scope for extension\. Richer feature representations could capture density structure not represented by low\-order moments, while alternative kernels, locally adaptive bandwidths, and multivariate density estimation provide natural extensions of the framework\. TheL1L^\{1\}approximation property of finite Gaussian mixtures further provides a principled basis for extending the training distribution over increasingly broad classes of density\-estimation tasks without explicitly enumerating individual distribution families\. Of particular practical interest are finite ensemble outputs, for which the underlying distributional shape and the appropriate amount of smoothing can vary from one forecast instance to another\. In such settings, a reusable selector that adapts its bandwidth to each observed ensemble while learning from a broad population of density\-estimation tasks provides a natural alternative to bandwidth rules whose functional form is fixed in advance\.

## 5Acknowledgments

This work was supported by the National Natural Science Foundation of China \(NSFC; Grant No\. 42450196\)\.

## Appendix

## Appendix AClassical Bandwidth Selectors

This appendix describes the three classical bandwidth selectors used throughout the experimental comparisons: Silverman’s rule, the Sheather–Jones selector, and least\-squares cross\-validation\. The same implementations are used in all experiments\. In the bounded\-support experiments, each selector is applied to the observed sample on its original scale, and the resulting bandwidth is subsequently used in the truncated\-and\-renormalized KDE defined in Section[2\.3](https://arxiv.org/html/2608.20445#S2.SS3)\. The corresponding location and scale transformation properties are established separately in Appendix[B](https://arxiv.org/html/2608.20445#A2)\.

### A\.1Silverman’s Rule

Silverman’s rule is a normal\-reference bandwidth selector that combines the usual sample standard deviation with a robust scale estimate based on the interquartile range\([18](https://arxiv.org/html/2608.20445#bib.bib1)\)\. For an observed sample

S=\(x1,…,xn\),S=\(x\_\{1\},\\ldots,x\_\{n\}\),\(A\.1\)the bandwidth used in the experiments is

hS\(S\)=0\.9min\{sd\(S\),IQR⁡\(S\)1\.34\}n−1/5,h\_\{\\mathrm\{S\}\}\(S\)=0\.9\\min\\left\\\{\\operatorname\{sd\}\(S\),\\frac\{\\operatorname\{IQR\}\(S\)\}\{1\.34\}\\right\\\}n^\{\-1/5\},\(A\.2\)where

sd⁡\(S\)=\[1n−1​∑i=1n\(xi−x¯\)2\]1/2\\operatorname\{sd\}\(S\)=\\left\[\\frac\{1\}\{n\-1\}\\sum\_\{i=1\}^\{n\}\(x\_\{i\}\-\\bar\{x\}\)^\{2\}\\right\]^\{1/2\}\(A\.3\)and

IQR⁡\(S\)=Q0\.75​\(S\)−Q0\.25​\(S\)\.\\operatorname\{IQR\}\(S\)=Q\_\{0\.75\}\(S\)\-Q\_\{0\.25\}\(S\)\.\(A\.4\)The minimum of the two scale estimates limits the influence of extreme observations while retaining the familiarn−1/5n^\{\-1/5\}dependence of the normal\-reference bandwidth\. The method is available in closed form and requires no sample\-specific numerical optimization\.

### A\.2Sheather–Jones Selector

The Sheather–Jones method is a data\-based plug\-in selector derived from the asymptotic mean integrated squared error of the KDE\([17](https://arxiv.org/html/2608.20445#bib.bib6)\)\. For a kernelκ\\kappa, define

R⁡\(κ\)=∫ℝκ​\(u\)2​𝑑u,μ2​\(κ\)=∫ℝu2​κ​\(u\)​𝑑u,R\(\\kappa\)=\\int\_\{\\mathbb\{R\}\}\\kappa\(u\)^\{2\}\\,du,\\qquad\\mu\_\{2\}\(\\kappa\)=\\int\_\{\\mathbb\{R\}\}u^\{2\}\\kappa\(u\)\\,du,\(A\.5\)and let

R⁡\(f′′\)=∫ℝ\{f′′​\(x\)\}2​𝑑xR\(f^\{\\prime\\prime\}\)=\\int\_\{\\mathbb\{R\}\}\\\{f^\{\\prime\\prime\}\(x\)\\\}^\{2\}\\,dx\(A\.6\)denote the roughness of the unknown density\. The leading asymptotic criterion has the form

AMISE⁡\(h\)=R⁡\(κ\)n​h\+h44​μ2​\(κ\)2​R​\(f′′\)\.\\operatorname\{AMISE\}\(h\)=\\frac\{R\(\\kappa\)\}\{nh\}\+\\frac\{h^\{4\}\}\{4\}\\mu\_\{2\}\(\\kappa\)^\{2\}R\(f^\{\\prime\\prime\}\)\.\(A\.7\)IfR⁡\(f′′\)R\(f^\{\\prime\\prime\}\)were known, minimizing this expression would give

hAMISE=\[R⁡\(κ\)n​μ2​\(κ\)2​R​\(f′′\)\]1/5\.h\_\{\\operatorname\{AMISE\}\}=\\left\[\\frac\{R\(\\kappa\)\}\{n\\mu\_\{2\}\(\\kappa\)^\{2\}R\(f^\{\\prime\\prime\}\)\}\\right\]^\{1/5\}\.\(A\.8\)For the standard Gaussian kernel,

R⁡\(κ\)=12​π,μ2​\(κ\)=1\.R\(\\kappa\)=\\frac\{1\}\{2\\sqrt\{\\pi\}\},\\qquad\\mu\_\{2\}\(\\kappa\)=1\.\(A\.9\)
The quantityR⁡\(f′′\)R\(f^\{\\prime\\prime\}\)is unknown and must be estimated from the sample\. The Sheather–Jones construction uses pilot estimates of density\-derivative functionals and determines the final bandwidth through a multistage plug\-in procedure\. In its solve\-the\-equation form, the defining relation can be represented schematically as

h=\[R⁡\(κ\)n​μ2​\(κ\)2​R^p⁡\(h\)​\(f′′\)\]1/5,h\>0,h=\\left\[\\frac\{R\(\\kappa\)\}\{n\\mu\_\{2\}\(\\kappa\)^\{2\}\\widehat\{R\}\_\{p\(h\)\}\(f^\{\\prime\\prime\}\)\}\\right\]^\{1/5\},\\qquad h\>0,\(A\.10\)wherep⁡\(h\)p\(h\)is a pilot bandwidth andR^p⁡\(h\)​\(f′′\)\\widehat\{R\}\_\{p\(h\)\}\(f^\{\\prime\\prime\}\)is the associated estimate of the roughness functional\. Because the pilot estimate depends on the final bandwidth throughp⁡\(h\)p\(h\), equation \([A\.10](https://arxiv.org/html/2608.20445#A1.E10)\) must be solved numerically\.

In the experiments, the Sheather–Jones bandwidth is computed using

bw\.SJ\(S, method = "ste"\)\(A\.11\)inR, where"ste"denotes the solve\-the\-equation version\. Its selected positive solution is denoted by

hSJ​\(S\)\.h\_\{\\mathrm\{SJ\}\}\(S\)\.\(A\.12\)The implementation uses pilot estimation of density derivatives and numerical root finding\. As with any numerical implementation, the reported value may be affected by floating\-point arithmetic, discretization, and root\-finding tolerance\.

### A\.3Least\-Squares Cross\-Validation

Least\-squares cross\-validation selects the bandwidth by minimizing a sample\-based estimate of integrated squared error\([15](https://arxiv.org/html/2608.20445#bib.bib3);[2](https://arxiv.org/html/2608.20445#bib.bib5)\)\. For a candidate bandwidthh\>0h\>0, let

qh​\(x,S\)=1n​h​∑i=1nϕ⁡\(x−xih\)q\_\{h\}\(x;S\)=\\frac\{1\}\{nh\}\\sum\_\{i=1\}^\{n\}\\phi\\left\(\\frac\{x\-x\_\{i\}\}\{h\}\\right\)\(A\.13\)be the Gaussian KDE, and define the leave\-one\-out estimator

qh,−i​\(x,S\)=1\(n−1\)​h​∑j=1j≠inϕ⁡\(x−xjh\)\.q\_\{h,\-i\}\(x;S\)=\\frac\{1\}\{\(n\-1\)h\}\\sum\_\{\\begin\{subarray\}\{c\}j=1\\\\ j\\neq i\\end\{subarray\}\}^\{n\}\\phi\\left\(\\frac\{x\-x\_\{j\}\}\{h\}\\right\)\.\(A\.14\)The least\-squares cross\-validation criterion is

CV⁡\(h;S\)=∫ℝqh​\(x,S\)2​𝑑x−2n​∑i=1nqh,−i​\(xi,S\)\.\\operatorname\{CV\}\(h;S\)=\\int\_\{\\mathbb\{R\}\}q\_\{h\}\(x;S\)^\{2\}\\,dx\-\\frac\{2\}\{n\}\\sum\_\{i=1\}^\{n\}q\_\{h,\-i\}\(x\_\{i\};S\)\.\(A\.15\)The term∫f​\(x\)2​𝑑x\\int f\(x\)^\{2\}\\,dxin the integrated squared error does not depend onhhand is therefore omitted from the criterion\.

For the Gaussian kernel, both terms in \([A\.15](https://arxiv.org/html/2608.20445#A1.E15)\) can be evaluated directly from pairwise sample differences\. Let

ϕ2​\(u\)=14​π​exp⁡\(−u24\)\\phi\_\{\\sqrt\{2\}\}\(u\)=\\frac\{1\}\{\\sqrt\{4\\pi\}\}\\exp\\left\(\-\\frac\{u^\{2\}\}\{4\}\\right\)\(A\.16\)denote the density of a centered Gaussian variable with variance22\. Then

CV⁡\(h,S\)\\displaystyle\\operatorname\{CV\}\(h;S\)=1n2​h​∑i=1n∑j=1nϕ2​\(xi−xjh\)\\displaystyle=\\frac\{1\}\{n^\{2\}h\}\\sum\_\{i=1\}^\{n\}\\sum\_\{j=1\}^\{n\}\\phi\_\{\\sqrt\{2\}\}\\left\(\\frac\{x\_\{i\}\-x\_\{j\}\}\{h\}\\right\)\(A\.17\)−2n⁡\(n−1\)​h∑i=1n∑j=1j≠inϕ\(xi−xjh\)\.\\displaystyle\-\\frac\{2\}\{n\(n\-1\)h\}\\sum\_\{i=1\}^\{n\}\\sum\_\{\\begin\{subarray\}\{c\}j=1\\\\ j\\neq i\\end\{subarray\}\}^\{n\}\\phi\\left\(\\frac\{x\_\{i\}\-x\_\{j\}\}\{h\}\\right\)\.
Rather than minimizing the criterion over allh\>0h\>0, the experimental implementation evaluates it over6060candidate bandwidths logarithmically spaced relative to Silverman’s bandwidth\. The candidate grid is

𝒢⁡\(S\)=\{h0​\(S\),…,h59​\(S\)\},\\mathcal\{G\}\(S\)=\\left\\\{h\_\{0\}\(S\),\\ldots,h\_\{59\}\(S\)\\right\\\},\(A\.18\)where

hk\(S\)=0\.2hS\(S\)\(50\.2\)k/59,k=0,…,59\.h\_\{k\}\(S\)=0\.2h\_\{\\mathrm\{S\}\}\(S\)\\left\(\\frac\{5\}\{0\.2\}\\right\)^\{k/59\},\\qquad k=0,\\ldots,59\.\(A\.19\)Thus,

h0​\(S\)=0\.2​hS​\(S\),h59​\(S\)=5​hS​\(S\),h\_\{0\}\(S\)=0\.2h\_\{\\mathrm\{S\}\}\(S\),\\qquad h\_\{59\}\(S\)=5h\_\{\\mathrm\{S\}\}\(S\),\(A\.20\)and adjacent grid points satisfy

hk\+1​\(S\)hk​\(S\)=\(50\.2\)1/59\.\\frac\{h\_\{k\+1\}\(S\)\}\{h\_\{k\}\(S\)\}=\\left\(\\frac\{5\}\{0\.2\}\\right\)^\{1/59\}\.\(A\.21\)The grid therefore has constant multiplicative spacing rather than constant additive spacing, providing comparable relative resolution across small and large bandwidths\.

The selected LSCV bandwidth is

hLSCV​\(S\)∈arg​minh∈𝒢⁡\(S\)​CV​\(h,S\)\.h\_\{\\mathrm\{LSCV\}\}\(S\)\\in\\operatorname\*\{arg\\,min\}\_\{h\\in\\mathcal\{G\}\(S\)\}\\operatorname\{CV\}\(h;S\)\.\(A\.22\)Defining the search interval relative tohS​\(S\)h\_\{\\mathrm\{S\}\}\(S\)also ensures that the complete candidate grid rescales with the observations, which is used in the scale\-equivariance argument in Appendix[B](https://arxiv.org/html/2608.20445#A2)\.

### A\.4Use in the Bounded\-Support Experiments

The definitions above determine the bandwidths before bounded\-support evaluation\. For a sample observed on\[A,B\]\[A,B\], each classical selector is applied directly to that sample on its original scale, producing

hb​\(S\),b∈\{S,SJ,LSCV\}\.h\_\{b\}\(S\),\\qquad b\\in\\\{\\mathrm\{S\},\\mathrm\{SJ\},\\mathrm\{LSCV\}\\\}\.\(A\.23\)The resulting bandwidth is then inserted into the truncated\-and\-renormalized estimator

q\[A,B\],hb​\(S\)​\(x,S\)=qhb​\(S\)​\(x,S\)​𝟏\[A,B\]​\(x\)∫ABqhb​\(S\)​\(u,S\)​𝑑u\.q\_\{\[A,B\],h\_\{b\}\(S\)\}\(x;S\)=\\frac\{q\_\{h\_\{b\}\(S\)\}\(x;S\)\\mathbf\{1\}\_\{\[A,B\]\}\(x\)\}\{\\displaystyle\\int\_\{A\}^\{B\}q\_\{h\_\{b\}\(S\)\}\(u;S\)\\,du\}\.\(A\.24\)Accordingly, all methods are evaluated using the same bounded\-support density representation and the same logarithmic score criterion\. The bounded\-support transformation changes the density used for evaluation, but it does not alter the bandwidth\-selection algorithms specified above\.

## Appendix BScale Equivariance of the Gaussian\-Comparison Bandwidth Selectors

This appendix establishes the location and scale transformation properties used in the Gaussian comparison\. Let

S=\(x1,…,xn\)S=\(x\_\{1\},\\ldots,x\_\{n\}\)\(B\.1\)be an observed sample, and define its affine transformation by

Sδ,c=δ\+c​S=\(δ\+c​x1,…,δ\+c​xn\),δ∈ℝ,c\>0\.S\_\{\\delta,c\}=\\delta\+cS=\(\\delta\+cx\_\{1\},\\ldots,\\delta\+cx\_\{n\}\),\\qquad\\delta\\in\\mathbb\{R\},\\quad c\>0\.\(B\.2\)A bandwidth selectorhbh\_\{b\}is translation invariant and scale equivariant if

hb​\(Sδ,c\)=c​hb​\(S\)\.h\_\{b\}\(S\_\{\\delta,c\}\)=c\\,h\_\{b\}\(S\)\.\(B\.3\)Thus, the bandwidth itself is not scale invariant: it changes in direct proportion to the scale of the data\. Dimensionless quantities such ashb​\(S\)/sd⁡\(S\)h\_\{b\}\(S\)/\\operatorname\{sd\}\(S\)are scale invariant\.

The derivations below concern the exact defining constructions of the selectors\. Numerical implementations may deviate from these exact identities because of floating\-point arithmetic, numerical root finding, binning, or the treatment of tied minima\.

### B\.1Affine Transformation of the Gaussian KDE and Logarithmic Score

For a bandwidthh\>0h\>0, the Gaussian KDE based onSSis

qh​\(x,S\)=1n​h​∑i=1nϕ⁡\(x−xih\),q\_\{h\}\(x;S\)=\\frac\{1\}\{nh\}\\sum\_\{i=1\}^\{n\}\\phi\\left\(\\frac\{x\-x\_\{i\}\}\{h\}\\right\),\(B\.4\)whereϕ\\phiis the standard\-normal density\. Consider the transformed sampleSδ,cS\_\{\\delta,c\}and the transformed bandwidthc​hch\. For anyy∈ℝy\\in\\mathbb\{R\},

qc​h​\(δ\+c​y,Sδ,c\)\\displaystyle q\_\{ch\}\(\\delta\+cy;S\_\{\\delta,c\}\)=1n​c​h​∑i=1nϕ⁡\(δ\+c​y−\(δ\+c​xi\)c​h\)\\displaystyle=\\frac\{1\}\{nch\}\\sum\_\{i=1\}^\{n\}\\phi\\left\(\\frac\{\\delta\+cy\-\(\\delta\+cx\_\{i\}\)\}\{ch\}\\right\)\(B\.5\)=1n​c​h​∑i=1nϕ⁡\(y−xih\)\\displaystyle=\\frac\{1\}\{nch\}\\sum\_\{i=1\}^\{n\}\\phi\\left\(\\frac\{y\-x\_\{i\}\}\{h\}\\right\)=1c​qh​\(y,S\)\.\\displaystyle=\\frac\{1\}\{c\}q\_\{h\}\(y;S\)\.Hence,

qc​h​\(δ\+c​y,Sδ,c\)=c−1​qh​\(y,S\)\.q\_\{ch\}\(\\delta\+cy;S\_\{\\delta,c\}\)=c^\{\-1\}q\_\{h\}\(y;S\)\.\(B\.6\)
Letffbe a density onℝ\\mathbb\{R\}and let

fδ,c​\(x\)=1c​f​\(x−δc\)f\_\{\\delta,c\}\(x\)=\\frac\{1\}\{c\}f\\left\(\\frac\{x\-\\delta\}\{c\}\\right\)\(B\.7\)be the density obtained by applying the transformationx↦δ\+c​xx\\mapsto\\delta\+cx\. The logarithmic score criterion for the transformed problem is

ℒ\(ch;Sδ,c,fδ,c\)=−∫ℝlog2qc​h\(x;Sδ,c\)fδ,c\(x\)dx\.\\mathcal\{L\}\\bigl\(ch;S\_\{\\delta,c\},f\_\{\\delta,c\}\\bigr\)=\-\\int\_\{\\mathbb\{R\}\}\\log\_\{2\}q\_\{ch\}\(x;S\_\{\\delta,c\}\)f\_\{\\delta,c\}\(x\)\\,dx\.\(B\.8\)Using the change of variablesx=δ\+c​yx=\\delta\+cy, so thatd​x=c​d​ydx=c\\,dy, gives

ℒ⁡\(c​h,Sδ,c,fδ,c\)\\displaystyle\\mathcal\{L\}\\bigl\(ch;S\_\{\\delta,c\},f\_\{\\delta,c\}\\bigr\)=−∫ℝlog2qc​h\(δ\+cy;Sδ,c\)f\(y\)dy\\displaystyle=\-\\int\_\{\\mathbb\{R\}\}\\log\_\{2\}q\_\{ch\}\(\\delta\+cy;S\_\{\\delta,c\}\)f\(y\)\\,dy\(B\.9\)=−∫ℝlog2\[1cqh\(y;S\)\]f\(y\)dy\\displaystyle=\-\\int\_\{\\mathbb\{R\}\}\\log\_\{2\}\\left\[\\frac\{1\}\{c\}q\_\{h\}\(y;S\)\\right\]f\(y\)\\,dy=−∫ℝlog2qh\(y;S\)f\(y\)dy\+log2c\\displaystyle=\-\\int\_\{\\mathbb\{R\}\}\\log\_\{2\}q\_\{h\}\(y;S\)f\(y\)\\,dy\+\\log\_\{2\}c=ℒ⁡\(h,S,f\)\+log2⁡c\.\\displaystyle=\\mathcal\{L\}\(h;S,f\)\+\\log\_\{2\}c\.Therefore,

ℒ⁡\(c​h,Sδ,c,fδ,c\)=ℒ⁡\(h,S,f\)\+log2⁡c\.\\boxed\{\\mathcal\{L\}\\bigl\(ch;S\_\{\\delta,c\},f\_\{\\delta,c\}\\bigr\)=\\mathcal\{L\}\(h;S,f\)\+\\log\_\{2\}c\.\}\(B\.10\)
The additive termlog2⁡c\\log\_\{2\}cdoes not depend on the bandwidth\. Consequently, the task\-specific minimizing bandwidths transform as

arg​minh′\>0⁡ℒ​\(h′,Sδ,c,fδ,c\)=c​arg​minh\>0⁡ℒ​\(h,S,f\)\.\\operatorname\*\{arg\\,min\}\_\{h^\{\\prime\}\>0\}\\mathcal\{L\}\\bigl\(h^\{\\prime\};S\_\{\\delta,c\},f\_\{\\delta,c\}\\bigr\)=c\\operatorname\*\{arg\\,min\}\_\{h\>0\}\\mathcal\{L\}\(h;S,f\)\.\(B\.11\)
Moreover, if the two bandwidth selectorshb1h\_\{b\_\{1\}\}andhb2h\_\{b\_\{2\}\}are scale equivariant, their logarithmic score difference is invariant under the same affine transformation:

ℒ⁡\(hb1​\(Sδ,c\),Sδ,c,fδ,c\)−ℒ⁡\(hb2​\(Sδ,c\),Sδ,c,fδ,c\)\\displaystyle\\mathcal\{L\}\\bigl\(h\_\{b\_\{1\}\}\(S\_\{\\delta,c\}\);S\_\{\\delta,c\},f\_\{\\delta,c\}\\bigr\)\-\\mathcal\{L\}\\bigl\(h\_\{b\_\{2\}\}\(S\_\{\\delta,c\}\);S\_\{\\delta,c\},f\_\{\\delta,c\}\\bigr\)\(B\.12\)=ℒ⁡\(hb1​\(S\),S,f\)−ℒ⁡\(hb2​\(S\),S,f\)\.\\displaystyle=\\mathcal\{L\}\\bigl\(h\_\{b\_\{1\}\}\(S\);S,f\\bigr\)\-\\mathcal\{L\}\\bigl\(h\_\{b\_\{2\}\}\(S\);S,f\\bigr\)\.

### B\.2Gaussian Amortized Selector

In the Gaussian comparison, the amortized selector is parameterized as

hθ​\(S\)=sd⁡\(S\)​Fθ​\(n\),h\_\{\\theta\}\(S\)=\\operatorname\{sd\}\(S\)F\_\{\\theta\}\(n\),\(B\.13\)whereFθ​\(n\)\>0F\_\{\\theta\}\(n\)\>0is dimensionless and depends only on the sample size\.

The sample mean ofSδ,cS\_\{\\delta,c\}is

x¯δ,c=δ\+c​x¯\.\\overline\{x\}\_\{\\delta,c\}=\\delta\+c\\overline\{x\}\.\(B\.14\)Therefore,

sd⁡\(Sδ,c\)\\displaystyle\\operatorname\{sd\}\(S\_\{\\delta,c\}\)=\[1n−1​∑i=1n\{\(δ\+c​xi\)−\(δ\+c​x¯\)\}2\]1/2\\displaystyle=\\left\[\\frac\{1\}\{n\-1\}\\sum\_\{i=1\}^\{n\}\\left\\\{\(\\delta\+cx\_\{i\}\)\-\(\\delta\+c\\overline\{x\}\)\\right\\\}^\{2\}\\right\]^\{1/2\}\(B\.15\)=\[c2n−1​∑i=1n\(xi−x¯\)2\]1/2\\displaystyle=\\left\[\\frac\{c^\{2\}\}\{n\-1\}\\sum\_\{i=1\}^\{n\}\(x\_\{i\}\-\\overline\{x\}\)^\{2\}\\right\]^\{1/2\}=c​sd⁡\(S\)\.\\displaystyle=c\\,\\operatorname\{sd\}\(S\)\.Since the sample size is unchanged,

hθ​\(Sδ,c\)\\displaystyle h\_\{\\theta\}\(S\_\{\\delta,c\}\)=sd⁡\(Sδ,c\)​Fθ​\(n\)\\displaystyle=\\operatorname\{sd\}\(S\_\{\\delta,c\}\)F\_\{\\theta\}\(n\)\(B\.16\)=c​sd⁡\(S\)​Fθ​\(n\)\\displaystyle=c\\,\\operatorname\{sd\}\(S\)F\_\{\\theta\}\(n\)=c​hθ​\(S\)\.\\displaystyle=c\\,h\_\{\\theta\}\(S\)\.Hence,

hθ​\(δ\+c​S\)=c​hθ​\(S\)\.\\boxed\{h\_\{\\theta\}\(\\delta\+cS\)=c\\,h\_\{\\theta\}\(S\)\.\}\(B\.17\)

### B\.3Silverman’s Rule

Silverman’s bandwidth is

hS\(S\)=0\.9min\{sd\(S\),IQR⁡\(S\)1\.34\}n−1/5\.h\_\{\\mathrm\{S\}\}\(S\)=0\.9\\min\\left\\\{\\operatorname\{sd\}\(S\),\\frac\{\\operatorname\{IQR\}\(S\)\}\{1\.34\}\\right\\\}n^\{\-1/5\}\.\(B\.18\)
As shown above,

sd⁡\(Sδ,c\)=c​sd⁡\(S\)\.\\operatorname\{sd\}\(S\_\{\\delta,c\}\)=c\\,\\operatorname\{sd\}\(S\)\.\(B\.19\)Forc\>0c\>0, empirical quantiles transform according to

Qp​\(Sδ,c\)=δ\+c​Qp​\(S\)\.Q\_\{p\}\(S\_\{\\delta,c\}\)=\\delta\+cQ\_\{p\}\(S\)\.\(B\.20\)It follows that

IQR⁡\(Sδ,c\)\\displaystyle\\operatorname\{IQR\}\(S\_\{\\delta,c\}\)=Q0\.75​\(Sδ,c\)−Q0\.25​\(Sδ,c\)\\displaystyle=Q\_\{0\.75\}\(S\_\{\\delta,c\}\)\-Q\_\{0\.25\}\(S\_\{\\delta,c\}\)\(B\.21\)=\[δ\+c​Q0\.75​\(S\)\]−\[δ\+c​Q0\.25​\(S\)\]\\displaystyle=\\bigl\[\\delta\+cQ\_\{0\.75\}\(S\)\\bigr\]\-\\bigl\[\\delta\+cQ\_\{0\.25\}\(S\)\\bigr\]=c​IQR⁡\(S\)\.\\displaystyle=c\\,\\operatorname\{IQR\}\(S\)\.
Substitution into Silverman’s formula gives

hS​\(Sδ,c\)\\displaystyle h\_\{\\mathrm\{S\}\}\(S\_\{\\delta,c\}\)=0\.9min\{csd\(S\),c​IQR⁡\(S\)1\.34\}n−1/5\\displaystyle=0\.9\\min\\left\\\{c\\,\\operatorname\{sd\}\(S\),\\frac\{c\\,\\operatorname\{IQR\}\(S\)\}\{1\.34\}\\right\\\}n^\{\-1/5\}\(B\.22\)=c0\.9min\{sd\(S\),IQR⁡\(S\)1\.34\}n−1/5\\displaystyle=c\\,0\.9\\min\\left\\\{\\operatorname\{sd\}\(S\),\\frac\{\\operatorname\{IQR\}\(S\)\}\{1\.34\}\\right\\\}n^\{\-1/5\}=c​hS​\(S\)\.\\displaystyle=c\\,h\_\{\\mathrm\{S\}\}\(S\)\.Therefore,

hS​\(δ\+c​S\)=c​hS​\(S\)\.\\boxed\{h\_\{\\mathrm\{S\}\}\(\\delta\+cS\)=c\\,h\_\{\\mathrm\{S\}\}\(S\)\.\}\(B\.23\)

### B\.4Sheather–Jones Selector

The scale equivariance of the Sheather–Jones selector follows from the homogeneity of its pilot derivative\-functional estimates\. The solve\-the\-equation construction can be expressed in terms of estimates of the form

ψ^r​\(p,S\)=1n2​pr\+1​∑i=1n∑j=1nL\(r\)​\(xi−xjp\),\\widehat\{\\psi\}\_\{r\}\(p;S\)=\\frac\{1\}\{n^\{2\}p^\{r\+1\}\}\\sum\_\{i=1\}^\{n\}\\sum\_\{j=1\}^\{n\}L^\{\(r\)\}\\left\(\\frac\{x\_\{i\}\-x\_\{j\}\}\{p\}\\right\),\(B\.24\)whereL\(r\)L^\{\(r\)\}is therrth derivative of a pilot kernel andp\>0p\>0is a pilot bandwidth\. Alternative normalizing conventions, such as usingn⁡\(n−1\)n\(n\-1\)and excluding diagonal terms, do not affect the scaling argument\.

Because pairwise differences are unaffected by translation and scale linearly withcc,

\(δ\+c​xi\)−\(δ\+c​xj\)c​p=xi−xjp\.\\frac\{\(\\delta\+cx\_\{i\}\)\-\(\\delta\+cx\_\{j\}\)\}\{cp\}=\\frac\{x\_\{i\}\-x\_\{j\}\}\{p\}\.\(B\.25\)Consequently,

ψ^r​\(c​p,Sδ,c\)\\displaystyle\\widehat\{\\psi\}\_\{r\}\(cp;S\_\{\\delta,c\}\)=1n2​\(c​p\)r\+1​∑i=1n∑j=1nL\(r\)​\(xi−xjp\)\\displaystyle=\\frac\{1\}\{n^\{2\}\(cp\)^\{r\+1\}\}\\sum\_\{i=1\}^\{n\}\\sum\_\{j=1\}^\{n\}L^\{\(r\)\}\\left\(\\frac\{x\_\{i\}\-x\_\{j\}\}\{p\}\\right\)\(B\.26\)=c−\(r\+1\)​ψ^r​\(p,S\)\.\\displaystyle=c^\{\-\(r\+1\)\}\\widehat\{\\psi\}\_\{r\}\(p;S\)\.Thus,

ψ^r​\(c​p,Sδ,c\)=c−\(r\+1\)​ψ^r​\(p,S\)\.\\boxed\{\\widehat\{\\psi\}\_\{r\}\(cp;S\_\{\\delta,c\}\)=c^\{\-\(r\+1\)\}\\widehat\{\\psi\}\_\{r\}\(p;S\)\.\}\(B\.27\)In particular,

ψ^4​\(c​p,Sδ,c\)=c−5​ψ^4​\(p,S\),\\widehat\{\\psi\}\_\{4\}\(cp;S\_\{\\delta,c\}\)=c^\{\-5\}\\widehat\{\\psi\}\_\{4\}\(p;S\),\(B\.28\)and

ψ^6​\(c​p,Sδ,c\)=c−7​ψ^6​\(p,S\)\.\\widehat\{\\psi\}\_\{6\}\(cp;S\_\{\\delta,c\}\)=c^\{\-7\}\\widehat\{\\psi\}\_\{6\}\(p;S\)\.\(B\.29\)
The solve\-the\-equation construction begins with a scale estimates⁡\(S\)s\(S\)satisfying

s⁡\(Sδ,c\)=c​s​\(S\)\.s\(S\_\{\\delta,c\}\)=c\\,s\(S\)\.\(B\.30\)For example, a scale estimate based on the minimum of the sample standard deviation and a rescaled interquartile range has this property\.

Let the initial pilot bandwidths be

p4\(S\)=C4s\(S\)n−1/7,p6\(S\)=C6s\(S\)n−1/9,p\_\{4\}\(S\)=C\_\{4\}\\,s\(S\)n^\{\-1/7\},\\qquad p\_\{6\}\(S\)=C\_\{6\}\\,s\(S\)n^\{\-1/9\},\(B\.31\)whereC4C\_\{4\}andC6C\_\{6\}are dimensionless constants\. Equation \([B\.30](https://arxiv.org/html/2608.20445#A2.E30)\) implies

p4​\(Sδ,c\)=c​p4​\(S\),p6​\(Sδ,c\)=c​p6​\(S\)\.p\_\{4\}\(S\_\{\\delta,c\}\)=c\\,p\_\{4\}\(S\),\\qquad p\_\{6\}\(S\_\{\\delta,c\}\)=c\\,p\_\{6\}\(S\)\.\(B\.32\)
Define the third\-derivative roughness pilot estimate by

R^3​\(S\)=−ψ^6​\(p6​\(S\),S\)\.\\widehat\{R\}\_\{3\}\(S\)=\-\\widehat\{\\psi\}\_\{6\}\\bigl\(p\_\{6\}\(S\);S\\bigr\)\.\(B\.33\)Using \([B\.29](https://arxiv.org/html/2608.20445#A2.E29)\) and \([B\.32](https://arxiv.org/html/2608.20445#A2.E32)\),

R^3​\(Sδ,c\)=c−7​R^3​\(S\)\.\\widehat\{R\}\_\{3\}\(S\_\{\\delta,c\}\)=c^\{\-7\}\\widehat\{R\}\_\{3\}\(S\)\.\(B\.34\)
The solve\-the\-equation pilot factor has the form

ηSJ​\(S\)=Cη​\[ψ^4​\(p4​\(S\),S\)R^3​\(S\)\]1/7,\\eta\_\{\\mathrm\{SJ\}\}\(S\)=C\_\{\\eta\}\\left\[\\frac\{\\widehat\{\\psi\}\_\{4\}\\bigl\(p\_\{4\}\(S\);S\\bigr\)\}\{\\widehat\{R\}\_\{3\}\(S\)\}\\right\]^\{1/7\},\(B\.35\)whereCηC\_\{\\eta\}is dimensionless\. From \([B\.28](https://arxiv.org/html/2608.20445#A2.E28)\), \([B\.32](https://arxiv.org/html/2608.20445#A2.E32)\), and \([B\.34](https://arxiv.org/html/2608.20445#A2.E34)\),

ηSJ​\(Sδ,c\)\\displaystyle\\eta\_\{\\mathrm\{SJ\}\}\(S\_\{\\delta,c\}\)=Cη​\[c−5​ψ^4​\(p4​\(S\),S\)c−7​R^3​\(S\)\]1/7\\displaystyle=C\_\{\\eta\}\\left\[\\frac\{c^\{\-5\}\\widehat\{\\psi\}\_\{4\}\\bigl\(p\_\{4\}\(S\);S\\bigr\)\}\{c^\{\-7\}\\widehat\{R\}\_\{3\}\(S\)\}\\right\]^\{1/7\}\(B\.36\)=c2/7​ηSJ​\(S\)\.\\displaystyle=c^\{2/7\}\\eta\_\{\\mathrm\{SJ\}\}\(S\)\.
For a candidate final bandwidthhh, define the corresponding pilot bandwidth by

pS​\(h\)=ηSJ​\(S\)​h5/7\.p\_\{S\}\(h\)=\\eta\_\{\\mathrm\{SJ\}\}\(S\)h^\{5/7\}\.\(B\.37\)Ifhhis replaced byc​hchandSSbySδ,cS\_\{\\delta,c\}, then

pSδ,c​\(c​h\)\\displaystyle p\_\{S\_\{\\delta,c\}\}\(ch\)=ηSJ​\(Sδ,c\)​\(c​h\)5/7\\displaystyle=\\eta\_\{\\mathrm\{SJ\}\}\(S\_\{\\delta,c\}\)\(ch\)^\{5/7\}\(B\.38\)=c2/7​ηSJ​\(S\)​c5/7​h5/7\\displaystyle=c^\{2/7\}\\eta\_\{\\mathrm\{SJ\}\}\(S\)c^\{5/7\}h^\{5/7\}=c​pS​\(h\)\.\\displaystyle=c\\,p\_\{S\}\(h\)\.
The Sheather–Jones bandwidth is defined as a positive solution of an equation of the form

h=\[Dnψ^4​\(pS​\(h\),S\)\]1/5,h=\\left\[\\frac\{D\_\{n\}\}\{\\widehat\{\\psi\}\_\{4\}\\bigl\(p\_\{S\}\(h\);S\\bigr\)\}\\right\]^\{1/5\},\(B\.39\)whereDnD\_\{n\}depends on the sample size and kernel constants but not on the scale of the observations\.

Suppose thathhsatisfies \([B\.39](https://arxiv.org/html/2608.20445#A2.E39)\) for the sampleSS\. For the transformed sample, evaluate the right\-hand side atc​hch\. From \([B\.27](https://arxiv.org/html/2608.20445#A2.E27)\) and \([B\.38](https://arxiv.org/html/2608.20445#A2.E38)\),

\[Dnψ^4​\(pSδ,c​\(c​h\),Sδ,c\)\]1/5\\displaystyle\\left\[\\frac\{D\_\{n\}\}\{\\widehat\{\\psi\}\_\{4\}\\bigl\(p\_\{S\_\{\\delta,c\}\}\(ch\);S\_\{\\delta,c\}\\bigr\)\}\\right\]^\{1/5\}\(B\.40\)=\[Dnψ^4​\(c​pS​\(h\),Sδ,c\)\]1/5\\displaystyle=\\left\[\\frac\{D\_\{n\}\}\{\\widehat\{\\psi\}\_\{4\}\\bigl\(c\\,p\_\{S\}\(h\);S\_\{\\delta,c\}\\bigr\)\}\\right\]^\{1/5\}=\[Dnc−5​ψ^4​\(pS​\(h\),S\)\]1/5\\displaystyle=\\left\[\\frac\{D\_\{n\}\}\{c^\{\-5\}\\widehat\{\\psi\}\_\{4\}\\bigl\(p\_\{S\}\(h\);S\\bigr\)\}\\right\]^\{1/5\}=c​\[Dnψ^4​\(pS​\(h\),S\)\]1/5\\displaystyle=c\\left\[\\frac\{D\_\{n\}\}\{\\widehat\{\\psi\}\_\{4\}\\bigl\(p\_\{S\}\(h\);S\\bigr\)\}\\right\]^\{1/5\}=c​h\.\\displaystyle=ch\.Thus, wheneverhhsolves the Sheather–Jones equation forSS, the bandwidthc​hchsolves the corresponding equation forSδ,cS\_\{\\delta,c\}\. Consequently, the sets of positive solutions satisfy

ℋSJ​\(Sδ,c\)=c​ℋSJ​\(S\)\.\\mathcal\{H\}\_\{\\mathrm\{SJ\}\}\(S\_\{\\delta,c\}\)=c\\,\\mathcal\{H\}\_\{\\mathrm\{SJ\}\}\(S\)\.\(B\.41\)When the defining equation has a unique positive solution, this reduces to

hSJ​\(δ\+c​S\)=c​hSJ​\(S\)\.\\boxed\{h\_\{\\mathrm\{SJ\}\}\(\\delta\+cS\)=c\\,h\_\{\\mathrm\{SJ\}\}\(S\)\.\}\(B\.42\)

### B\.5Least\-Squares Cross\-Validation on a Relative Grid

For a bandwidthhh, the least\-squares cross\-validation criterion is

CV⁡\(h;S\)=∫ℝqh​\(x,S\)2​𝑑x−2n​∑i=1nqh,−i​\(xi,S\),\\operatorname\{CV\}\(h;S\)=\\int\_\{\\mathbb\{R\}\}q\_\{h\}\(x;S\)^\{2\}\\,dx\-\\frac\{2\}\{n\}\\sum\_\{i=1\}^\{n\}q\_\{h,\-i\}\(x\_\{i\};S\),\(B\.43\)where

qh,−i​\(x,S\)=1\(n−1\)​h​∑j≠iϕ⁡\(x−xjh\)q\_\{h,\-i\}\(x;S\)=\\frac\{1\}\{\(n\-1\)h\}\\sum\_\{j\\neq i\}\\phi\\left\(\\frac\{x\-x\_\{j\}\}\{h\}\\right\)\(B\.44\)is the leave\-one\-out KDE\.

The KDE transformation in \([B\.6](https://arxiv.org/html/2608.20445#A2.E6)\) gives

qc​h​\(δ\+c​y,Sδ,c\)=c−1​qh​\(y,S\)\.q\_\{ch\}\(\\delta\+cy;S\_\{\\delta,c\}\)=c^\{\-1\}q\_\{h\}\(y;S\)\.\(B\.45\)The leave\-one\-out estimator satisfies the analogous relation

qc​h,−i​\(δ\+c​xi,Sδ,c\)=c−1​qh,−i​\(xi,S\)\.q\_\{ch,\-i\}\(\\delta\+cx\_\{i\};S\_\{\\delta,c\}\)=c^\{\-1\}q\_\{h,\-i\}\(x\_\{i\};S\)\.\(B\.46\)
For the integrated\-square term, the change of variablesx=δ\+c​yx=\\delta\+cygives

∫ℝqc​h​\(x,Sδ,c\)2​𝑑x\\displaystyle\\int\_\{\\mathbb\{R\}\}q\_\{ch\}\(x;S\_\{\\delta,c\}\)^\{2\}\\,dx=∫ℝqc​h​\(δ\+c​y,Sδ,c\)2​c​𝑑y\\displaystyle=\\int\_\{\\mathbb\{R\}\}q\_\{ch\}\(\\delta\+cy;S\_\{\\delta,c\}\)^\{2\}c\\,dy\(B\.47\)=∫ℝc−2​qh​\(y,S\)2​c​𝑑y\\displaystyle=\\int\_\{\\mathbb\{R\}\}c^\{\-2\}q\_\{h\}\(y;S\)^\{2\}c\\,dy=c−1​∫ℝqh​\(y,S\)2​dy\.\\displaystyle=c^\{\-1\}\\int\_\{\\mathbb\{R\}\}q\_\{h\}\(y;S\)^\{2\}\\,dy\.The leave\-one\-out term satisfies

2n​∑i=1nqc​h,−i​\(δ\+c​xi,Sδ,c\)\\displaystyle\\frac\{2\}\{n\}\\sum\_\{i=1\}^\{n\}q\_\{ch,\-i\}\(\\delta\+cx\_\{i\};S\_\{\\delta,c\}\)=c−1​2n​∑i=1nqh,−i​\(xi,S\)\.\\displaystyle=c^\{\-1\}\\frac\{2\}\{n\}\\sum\_\{i=1\}^\{n\}q\_\{h,\-i\}\(x\_\{i\};S\)\.\(B\.48\)Combining \([B\.47](https://arxiv.org/html/2608.20445#A2.E47)\) and \([B\.48](https://arxiv.org/html/2608.20445#A2.E48)\) yields

CV⁡\(c​h;Sδ,c\)=c−1​CV⁡\(h;S\)\.\\boxed\{\\operatorname\{CV\}\\bigl\(ch;S\_\{\\delta,c\}\\bigr\)=c^\{\-1\}\\operatorname\{CV\}\(h;S\)\.\}\(B\.49\)
The multiplicative factorc−1c^\{\-1\}is positive and independent of the bandwidth\. It therefore preserves the ordering of candidate bandwidths\.

In the present implementation, LSCV is evaluated over the grid

hk\(S\)=0\.2hS\(S\)\(50\.2\)k/59,k=0,…,59\.h\_\{k\}\(S\)=0\.2h\_\{\\mathrm\{S\}\}\(S\)\\left\(\\frac\{5\}\{0\.2\}\\right\)^\{k/59\},\\qquad k=0,\\ldots,59\.\(B\.50\)Because Silverman’s bandwidth satisfies \([B\.23](https://arxiv.org/html/2608.20445#A2.E23)\),

hk​\(Sδ,c\)\\displaystyle h\_\{k\}\(S\_\{\\delta,c\}\)=0\.2​hS​\(Sδ,c\)​\(50\.2\)k/59\\displaystyle=0\.2h\_\{\\mathrm\{S\}\}\(S\_\{\\delta,c\}\)\\left\(\\frac\{5\}\{0\.2\}\\right\)^\{k/59\}\(B\.51\)=c​hk​\(S\)\.\\displaystyle=c\\,h\_\{k\}\(S\)\.Thus, the candidate sets satisfy

𝒢⁡\(Sδ,c\)=c​𝒢​\(S\),\\mathcal\{G\}\(S\_\{\\delta,c\}\)=c\\,\\mathcal\{G\}\(S\),\(B\.52\)where

𝒢⁡\(S\)=\{h0​\(S\),…,h59​\(S\)\}\.\\mathcal\{G\}\(S\)=\\\{h\_\{0\}\(S\),\\ldots,h\_\{59\}\(S\)\\\}\.\(B\.53\)
For every corresponding grid point,

CV⁡\(hk​\(Sδ,c\);Sδ,c\)=c−1​CV⁡\(hk​\(S\);S\)\.\\operatorname\{CV\}\\bigl\(h\_\{k\}\(S\_\{\\delta,c\}\);S\_\{\\delta,c\}\\bigr\)=c^\{\-1\}\\operatorname\{CV\}\\bigl\(h\_\{k\}\(S\);S\\bigr\)\.\(B\.54\)All6060criterion values are therefore multiplied by the same positive constant, so the minimizing grid index is unchanged\. It follows that

hLSCV​\(δ\+c​S\)=c​hLSCV​\(S\)\.\\boxed\{h\_\{\\mathrm\{LSCV\}\}\(\\delta\+cS\)=c\\,h\_\{\\mathrm\{LSCV\}\}\(S\)\.\}\(B\.55\)If several grid points attain the same minimum, the corresponding minimizer sets transform in the same way\. A deterministic tie\-breaking rule based on the grid index also preserves the result\.

### B\.6Consequence for the Gaussian Comparison

Let

S0∼𝒩​\(0,1\)nS\_\{0\}\\sim\\mathcal\{N\}\(0,1\)^\{n\}\(B\.56\)and define

Sδ,c=δ\+c​S0∼𝒩​\(δ,c2\)n\.S\_\{\\delta,c\}=\\delta\+cS\_\{0\}\\sim\\mathcal\{N\}\(\\delta,c^\{2\}\)^\{n\}\.\(B\.57\)For every selector considered in the Gaussian comparison,

hb​\(Sδ,c\)=c​hb​\(S0\),b∈\{θ,S,SJ,LSCV\}\.h\_\{b\}\(S\_\{\\delta,c\}\)=c\\,h\_\{b\}\(S\_\{0\}\),\\qquad b\\in\\\{\\theta,\\mathrm\{S\},\\mathrm\{SJ\},\\mathrm\{LSCV\}\\\}\.\(B\.58\)Equation \([B\.10](https://arxiv.org/html/2608.20445#A2.E10)\) then gives

ℒ⁡\(hb​\(Sδ,c\),Sδ,c,fδ,c\)=ℒ⁡\(hb​\(S0\),S0,ϕ\)\+log2⁡c\.\\mathcal\{L\}\\bigl\(h\_\{b\}\(S\_\{\\delta,c\}\);S\_\{\\delta,c\},f\_\{\\delta,c\}\\bigr\)=\\mathcal\{L\}\\bigl\(h\_\{b\}\(S\_\{0\}\);S\_\{0\},\\phi\\bigr\)\+\\log\_\{2\}c\.\(B\.59\)
The absolute logarithmic scores for𝒩⁡\(δ,c2\)\\mathcal\{N\}\(\\delta,c^\{2\}\)therefore differ from those for𝒩⁡\(0,1\)\\mathcal\{N\}\(0,1\)by the common additive termlog2⁡c\\log\_\{2\}c\. The ordering of the selectors and all pairwise score differences are unchanged\. Standard\-normal sampling is consequently sufficient for the comparison of these bandwidth selectors throughout the Gaussian location–scale family\.

## Appendix CNeural Network Architecture and Training Procedure

### C\.1Gaussian Benchmark

For the Gaussian benchmark, the amortized selector uses the scale\-equivariant parameterization

hθ​\(S\)=sd⁡\(S\)​Fθ​\(n\),Fθ​\(n\)\>0,h\_\{\\theta\}\(S\)=\\operatorname\{sd\}\(S\)\\,F\_\{\\theta\}\(n\),\\qquad F\_\{\\theta\}\(n\)\>0,whereFθ​\(n\)F\_\{\\theta\}\(n\)is predicted from the sample size alone\. The predictor is a fully connected multilayer perceptron with layer widths

1⟶4⟶4⟶1,1\\longrightarrow 4\\longrightarrow 4\\longrightarrow 1,with ReLU activations in the hidden layers and a softplus output transformation\.

Training tasks are generated online\. At each optimization step,n∼Unif⁡\{5,…,256\}n\\sim\\operatorname\{Unif\}\\\{5,\\ldots,256\\\}is drawn, and a batch of 256 independent samples from𝒩​\(0,1\)n\\mathcal\{N\}\(0,1\)^\{n\}is generated\. The expected logarithmic score for each task is evaluated using 256\-point Gauss–Hermite quadrature\.

The network is optimized using Adam with learning rate10−310^\{\-3\}and gradient clipping at11, for at most 10,000 optimization steps\. The training seed is fixed at 1234\. Model selection is based on a fixed validation set over

\{5,10,15,20,30,50,75,100,150,200,256\},\\\{5,10,15,20,30,50,75,100,150,200,256\\\},with 512 independent tasks at each sample size\. Validation is performed every 500 steps\. Early stopping is applied after 3000 steps and is triggered after three consecutive validation evaluations without an improvement of at least10−410^\{\-4\}bits\. The checkpoint used in the Gaussian comparison was retained at optimization step 3500\.

### C\.2Bounded Multi\-Family Benchmark

For the bounded multi\-family benchmark, the amortized selector uses the five\-dimensional feature representation in \(4\)\. The predictor is a fully connected network with layer widths

5⟶128⟶128⟶1,5\\longrightarrow 128\\longrightarrow 128\\longrightarrow 1,with ReLU hidden activations and a softplus output transformation to ensure a positive bandwidth\.

As described in Section 3\.2 and Table 1, 10,000 parameter settings are generated in advance for each of the ten distribution families\. For each training task, a family is sampled with equal probability, a parameter setting is drawn from the corresponding family\-specific pool, and a fresh sample is generated from the resulting bounded distribution\. The sample size is drawn independently for each task fromn∼Unif⁡\{5,…,256\}n\\sim\\operatorname\{Unif\}\\\{5,\\ldots,256\\\}\. Each training batch contains 256 tasks, and the expected bounded logarithmic score for each task is approximated using 1024 independent Monte Carlo draws from the same target distribution\.

Training uses Adam with learning rate5×10−45\\times 10^\{\-4\}, gradient clipping at11, and an exponential moving average of the network parameters with decay0\.9990\.999\. Training is allowed to continue for at most 40,000 optimization steps, with training seed 2026\.

Model selection uses a fixed validation set containing 20 tasks for every integern∈\{5,…,256\}n\\in\\\{5,\\ldots,256\\\}, giving 5040 tasks in total\. The validation set is generated once with seed 777 and held fixed throughout training\. Validation is performed every 1000 steps\. Early stopping is not applied before 12,000 steps and is triggered after five consecutive validation evaluations without an improvement of at least5×10−55\\times 10^\{\-5\}bits\. Training reached the maximum of 40,000 optimization steps, at which the best checkpoint was retained\. The corresponding exponential\-moving\-average weights are used for subsequent evaluation\.

### C\.3GMM Benchmarks

TheK=8K=8,K=16K=16, andK=32K=32GMM\-trained selectors use the same feature representation, network architecture, bounded logarithmic score objective, and training procedure as the multi\-family selector\.

The main difference is the generation of training tasks\. Rather than sampling from a fixed parameter pool, a fresh GMM is independently generated for every training task using the corresponding generator defined in Section 3\.3\. The sample size is drawn independently for each task fromn∼Unif⁡\{5,…,256\}n\\sim\\operatorname\{Unif\}\\\{5,\\ldots,256\\\}, and both the observed sample and an independent Monte Carlo evaluation sample of size 1024 are generated from the same bounded GMM\.

Separate networks are trained forK=8K=8,K=16K=16, andK=32K=32\. The optimizer, batch size, exponential moving average, validation schedule, early\-stopping rule, and random seeds are the same as in the multi\-family experiment\. For each value ofKK, model selection is based on an independently generated fixed validation set\. All three training runs reached the maximum of 40,000 optimization steps without triggering early stopping\. The best checkpoints were retained at steps 39,000, 39,000, and 40,000 forK=8K=8,K=16K=16, andK=32K=32, respectively, and their exponential\-moving\-average weights are used in the final comparisons\.

### C\.4Evaluation and Bootstrap Procedure

All bandwidth\-comparison experiments are evaluated on the sample\-size grid

\{5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,25,30,35,40,50,60,75,90,110,130,150,175,200,225,256\}\.\\begin\{split\}\\\{&5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,\\\\ &25,30,35,40,50,60,75,90,110,130,150,175,200,225,256\\\}\.\\end\{split\}
For the Gaussian benchmark, 30,000 independent test samples are generated at each sample size\. For the bounded multi\-family benchmark, 3000 independent test tasks are generated from each of the ten families at each sample size\. For the bounded GMM benchmark, 30,000 independently generatedK=32K=32GMM tasks are used at each sample size\. The cross\-distribution comparison on the multi\-family benchmark uses the same balanced design of 3000 tasks per family\.

For the bounded bandwidth\-comparison experiments, each task uses 1024 independent Monte Carlo evaluation draws\. Within each experiment, all bandwidth selectors are evaluated using the same task\-specific samples and evaluation draws\.

Uncertainty bands are obtained using 1000 paired nonparametric bootstrap replicates\. At each sample size, 30,000 task indices are resampled with replacement, and the same indices are used for all methods\. The lower and upper limits are the empirical 5th and 95th percentiles of the bootstrap means, giving 90% bootstrap intervals\.

## References

- J\. M\. BernardoExpected information as expected utility\.The Annals of Statistics7\(3\),pp\. 686–690\.Cited by:[§1](https://arxiv.org/html/2608.20445#S1.p3.1),[§2\.2](https://arxiv.org/html/2608.20445#S2.SS2.p2.1)\.
- Bowman \(1984\)A\. W\. BowmanAn alternative method of cross\-validation for the smoothing of density estimates\.Biometrika71\(2\),pp\. 353–360\.Cited by:[§A\.3](https://arxiv.org/html/2608.20445#A1.SS3.p1.1),[§1](https://arxiv.org/html/2608.20445#S1.p2.1)\.
- Bröcker and Smith \(2007\)J\. Bröcker and L\. A\. SmithScoring probabilistic forecasts: the importance of being proper\.Weather and Forecasting22\(2\),pp\. 382–388\.Cited by:[§1](https://arxiv.org/html/2608.20445#S1.p3.1),[§2\.2](https://arxiv.org/html/2608.20445#S2.SS2.p1.1),[§2\.2](https://arxiv.org/html/2608.20445#S2.SS2.p2.1)\.
- Du \(2021\)H\. DuBeyond strictly proper scoring rules: the importance of being local\.Weather and Forecasting36\(2\),pp\. 457–468\.Cited by:[§1](https://arxiv.org/html/2608.20445#S1.p3.1),[§2\.2](https://arxiv.org/html/2608.20445#S2.SS2.p1.1),[§2\.2](https://arxiv.org/html/2608.20445#S2.SS2.p2.1),[§3\.1](https://arxiv.org/html/2608.20445#S3.SS1.p4.1),[§3](https://arxiv.org/html/2608.20445#S3.p2.1)\.
- Gneitinget al\.\(2007\)T\. Gneiting, F\. Balabdaoui, and A\. E\. RafteryProbabilistic forecasts, calibration and sharpness\.Journal of the Royal Statistical Society: Series B \(Statistical Methodology\)69\(2\),pp\. 243–268\.External Links:[Document](https://dx.doi.org/10.1111/j.1467-9868.2007.00587.x)Cited by:[§2\.2](https://arxiv.org/html/2608.20445#S2.SS2.p1.1)\.
- Gneiting and Raftery \(2007\)T\. Gneiting and A\. E\. RafteryStrictly proper scoring rules, prediction, and estimation\.Journal of the American Statistical Association102\(477\),pp\. 359–378\.Cited by:[§2\.2](https://arxiv.org/html/2608.20445#S2.SS2.p1.1)\.
- Good \(1952\)I\. J\. GoodRational decisions\.Journal of the Royal Statistical Society: Series B \(Methodological\)14\(1\),pp\. 107–114\.Cited by:[§2\.2](https://arxiv.org/html/2608.20445#S2.SS2.p2.1)\.
- Green and Richardson \(2001\)P\. J\. Green and S\. RichardsonModelling heterogeneity with and without the Dirichlet process\.Scandinavian Journal of Statistics28\(2\),pp\. 355–375\.External Links:[Document](https://dx.doi.org/10.1111/1467-9469.00242)Cited by:[§3\.3](https://arxiv.org/html/2608.20445#S3.SS3.p3.2)\.
- Hall \(1983\)P\. HallLarge sample optimality of least squares cross\-validation in density estimation\.The Annals of Statistics11\(4\),pp\. 1156–1174\.Cited by:[§1](https://arxiv.org/html/2608.20445#S1.p2.1)\.
- Jewson and Rossell \(2022\)J\. Jewson and D\. RossellGeneral Bayesian loss function selection and the use of improper models\.Journal of the Royal Statistical Society: Series B \(Statistical Methodology\)84\(5\),pp\. 1640–1665\.External Links:[Document](https://dx.doi.org/10.1111/rssb.12553)Cited by:[§2\.2](https://arxiv.org/html/2608.20445#S2.SS2.p2.1)\.
- Malsiner\-Walliet al\.\(2016\)G\. Malsiner\-Walli, S\. Frühwirth\-Schnatter, and B\. GrünModel\-based clustering based on sparse finite Gaussian mixtures\.Statistics and Computing26\(1–2\),pp\. 303–324\.External Links:[Document](https://dx.doi.org/10.1007/s11222-014-9500-2)Cited by:[§3\.3](https://arxiv.org/html/2608.20445#S3.SS3.p5.1)\.
- Nguyenet al\.\(2023\)T\. Nguyen, F\. Chamroukhi, H\. D\. Nguyen, and G\. J\. McLachlanApproximation of probability density functions via location\-scale finite mixtures in Lebesgue spaces\.Communications in Statistics – Theory and Methods52\(14\),pp\. 5048–5059\.External Links:[Document](https://dx.doi.org/10.1080/03610926.2021.2002360)Cited by:[§3\.3](https://arxiv.org/html/2608.20445#S3.SS3.p2.2)\.
- Richardson and Green \(1997\)S\. Richardson and P\. J\. GreenOn Bayesian analysis of mixtures with an unknown number of components\.Journal of the Royal Statistical Society: Series B \(Methodological\)59\(4\),pp\. 731–792\.External Links:[Document](https://dx.doi.org/10.1111/1467-9868.00095)Cited by:[§3\.3](https://arxiv.org/html/2608.20445#S3.SS3.p5.1)\.
- Roulston and Smith \(2002\)M\. S\. Roulston and L\. A\. SmithEvaluating probabilistic forecasts using information theory\.Monthly Weather Review130\(6\),pp\. 1653–1660\.Cited by:[§2\.2](https://arxiv.org/html/2608.20445#S2.SS2.p2.1)\.
- Rudemo \(1982\)M\. RudemoEmpirical choice of histograms and kernel density estimators\.Scandinavian Journal of Statistics9\(2\),pp\. 65–78\.Cited by:[§A\.3](https://arxiv.org/html/2608.20445#A1.SS3.p1.1),[§1](https://arxiv.org/html/2608.20445#S1.p2.1)\.
- Scott \(2015\)D\. W\. ScottMultivariate density estimation: theory, practice, and visualization\.2nd edition,John Wiley & Sons,Hoboken, NJ\.Cited by:[§1](https://arxiv.org/html/2608.20445#S1.p1.1),[§1](https://arxiv.org/html/2608.20445#S1.p2.1),[§1](https://arxiv.org/html/2608.20445#S1.p3.1)\.
- Sheather and Jones \(1991\)S\. J\. Sheather and M\. C\. JonesA reliable data\-based bandwidth selection method for kernel density estimation\.Journal of the Royal Statistical Society: Series B \(Methodological\)53\(3\),pp\. 683–690\.Cited by:[§A\.2](https://arxiv.org/html/2608.20445#A1.SS2.p1.1),[§1](https://arxiv.org/html/2608.20445#S1.p2.1)\.
- Silverman \(1986\)B\. W\. SilvermanDensity estimation for statistics and data analysis\.Chapman and Hall,London\.Cited by:[§A\.1](https://arxiv.org/html/2608.20445#A1.SS1.p1.1),[§1](https://arxiv.org/html/2608.20445#S1.p1.1),[§1](https://arxiv.org/html/2608.20445#S1.p2.1),[§1](https://arxiv.org/html/2608.20445#S1.p3.1)\.
- Wand and Jones \(1995\)M\. P\. Wand and M\. C\. JonesKernel smoothing\.Chapman and Hall/CRC,London\.Cited by:[§1](https://arxiv.org/html/2608.20445#S1.p2.1),[§1](https://arxiv.org/html/2608.20445#S1.p3.1)\.