## 2 Materials and Methods
### Kernel Density Estimation  
[PERSON] [43] formulated univariate kernel density estimation (KDE) in its modern form including the smoothing parameter, that is, bandwidth \(h\), as:  
\[p_{h}(x)=\frac{1}{nh}\sum_{j=1}^{n}K\bigg{(}\frac{x-x_{j}}{h}\bigg{)}, \tag{4}\]where the kernel \(K\) is a non-negative function and \(x_{j}\) are the \(n\) data points. We use an analogously defined multivariate version of KDE [20; 44] with a bandwidth matrix \(\mathbf{S}\) as:  
\[p_{h}(\mathbf{x})=\frac{1}{n}\sum_{j=1}^{n}K_{h}(\mathbf{x}-\mathbf{x}_{j}), \tag{5}\]  
with the standard Gaussian kernel \(K_{h}(\mathbf{x})=(2\pi)^{-D/2}|\mathbf{S}|^{-1/2}e^{-\frac{1}{2}\mathbf{x}^{T}yellow \mathbf{S}^{-1}\mathbf{x}}\) and a diagonal bandwidth matrix as the covariance matrix \(\sqrt{S_{dd}}=n^{\frac{-1}{D+1}}h_{d}\) determined by [PERSON]'s rule [21]. Note that for small objects the estimator is smoothed more, due to an inverse relationship between \(n\) and \(\mathbf{S}_{dd}\). We refer to this estimate as Gaussian KDE (GKDE). GKDE is an effective, lightweight method of providing smoothed probability density estimates for point samples independently of discretization interval or data point count. However, GKDE provides no measure of uncertainty relative to its suggested point estimate, and hence, similarly to histograms, loses information about the relative reliability of different objects. Logistic Gaussian Process Density Estimation  
For objects with only a few pixels, it becomes important to explicitly quantify the uncertainty of the density estimate itself, which neither of the above methods can achieve. For instance, for the extreme case of just one pixel, the histogram becomes a delta distribution, and while GKDE provides a smoother estimate it still suggests this single noisy pixel alone to be highly informative of the content. Bayesian estimators, instead, have the ability to explicitly model uncertainty, and in the following we describe one practical alternative building on Gaussian Processes (GP). LGPDE, originally proposed by [PERSON] et al. [45], assigns a GP prior for the unnormalized logarithmic density \(f(\mathbf{x})\) so that \(\log p(\mathbf{x})=f(\mathbf{x})+C\) for any \(\mathbf{x}\), where \(C\) is a constant required for normalizing the density. The GP assigns a prior over the functions directly, so that for any finite collection of inputs their joint distribution is a multivariate normal, and conditioning on some pixel observations \(X\) we can then obtain the posterior distribution \(p(f|X)\) that captures the uncertainty of the estimator. Due to the logistic transform, there is no closed-form analytic expression for the posterior, but both Markov Chain Monte Carlo (MCMC) sampling [46] and Laplace approximation [22] can be used for inference. We will later evaluate both the Laplace approximation as well as an MCMC implementation using the No-U-Turn Hamiltonian Monte Carlo algorithm as provided in the Stan probabilistic programming environment [47]. We use the formulation of [PERSON] et al. [22] with explicit enumeration over discretized support axes for computing the normalization term \(C\). A prior term results from the logistic transform:  
\[\log p(\mathbf{f}|G,\theta)=\mathcal{N}(\mathbf{f}|H\mathbf{m},K+HMH^{T}), \tag{6}\]  
where \(\mathbf{f}\) is a latent function representing the density estimate surface being evaluated at points \(\mathbf{x}_{g}\) of the discretization grid \(G\), \(\theta\) denotes the hyperparameters of the prior and the GP kernel, and \(H(G)\) is a basis function that modulates the density to achieve finite support. For 2D densities we use the basis function \(H(\mathbf{x})=[x_{1},x_{1}^{2},x_{2},x_{2}^{2},x_{1}x_{2}]^{T}\). For a weakly informative prior, we parametrize a covariance adjustment of \(M=10^{2}I\) and a zero mean of \(\mathbf{m}=\mathbf{0}\). The kernel \(K=K(G)\) determines a covariance matrix based on a given covariance function and a chosen multivariate bin discretization expressed by \(G\).