fast & accurate gaussian kernel density estimation
for a set of input data points…
…we represent each point with a kernel function.
we then sum the kernels to form a continuous density estimate. here we focus on gaussian kernels which are commonly used in practice.
the spread of the kernels is determined by a bandwidth parameter for gaussian kernels this is just the standard deviation.
unfortunately direct calculation is expensive. for each of m measurement points we must sum the contributions of the n data points resulting in quadratic complexity.
so how might we achieve fast yet accurate approximate estimates? one approach is to discretize the problem by binning the data. kde of n data points then reduces to a signal processing task smoothing a grid of m bins.
with simple binning the weight of a data point is assigned to the nearest bin. if we perturb the points we can see how the weight is reassigned across bin boundaries
alternatively we can linearly interpolate weight across bins the weight is assigned proportionally providing some smoothing. we evaluate both of these binning methods.
now, a number of gaussian smoothing approaches exist for binned data.
we focus on linear time methods that first bin the data, an o(n) step, and then make a bounded number of passes over the binned grid, an o(m) step.
the first method we'll consider is box filtering which approximates gaussian smoothing by iteratively applying a uniform or box filter .
we sum and scale the input values in the filter window to produce an output bin value.
we apply this running sum across the input grid writing results into an output grid.
to perform another iteration we first swap the input and output grids and then apply the box filter again.
but note that the box filters can blur weight outside our original grid extent. box filter methods require grids with extra padding, which may increase memory use and running time.
we perform our second iteration…
…swap the grids again perform a third iteration…
…and arrive at our density estimate. box filters are fast and easy to implement but they do have important nuances.
we already saw how we must add padding to the grid, but in addition the filter width, which is a function of bandwidth, they can suffer from quantization error. here bandwidths of 0.1 and 0.15 both map to a box filter of width 3. the box filter results for the two will be indistinguishable.
similarly here we get a filter width of 5 bins for bandwidths of both 0.2 and 0.25
to address this, the extended box filter method adds fractional weight to the ends of the filter window .
as the bandwidth varies the extended box varies the filter response more smoothly.
finally let's consider an approximation developed in the early 90s by deriche , published in the signal processing literature. this method seems to have been overlooked by statistics and visualization researchers.
deriche derived a recursive filter approximation for the right half of a gaussian with parameters he determined via optimization.
we can reverse the equation to model the left half of a gaussian and then sum the two filter responses.
the method runs in linear time but does involve more arithmetic operations than box filters.
so now let's evaluate how these linear time methods fare in terms of accuracy.
we’ll start with an impulse test involving a single data point, corresponding to a single gaussian density.
we can compare direct calculation with box filters…
…and extended box filters. we see that both methods lead to underestimation of the peak and poor fit to the tails of the distribution.
meanwhile, the deriche method achieves high accuracy.
we must also inspect other bandwidth values…
…to more systematically measure the error.
for each kde method we measure the maximum pixel error for a 100 pixel tall density chart, plotting the results on a logarithmic scale. we do this across a range of bandwidths.
here are the error curves for simple binning on a 256 bin grid. we see oscillation in the box method due to filter quantization. the extended box method smooths this out. deriche's method outperforms the others and, at high bandwidth, it's over an order of magnitude better with a maximum error around one pixel.
switching to linear binning the deriche method improves with subpixel accuracy for most tested bandwidths
if we double the number of bins deriche improves further to near pixel perfect accuracy almost two orders of magnitude better than box filter methods.
ultimately we care about real world data, so here's an excerpt from the palmer penguins data set showing an estimated density of penguin body mass.
at bandwidth 200, all methods appear to perform reasonably well.
however at lower bandwidths we see some significant discrepancies. the box filter methods may misestimate peaks or erode local optima and in the language of algebraic vis these are hallucinators that is our choice of technique may mislead us with inaccurate visual features .
looking across bandwidths, here are the results for simple binning with 256 bins. again the deriche method provides the best performance.
accuracy improves with linear binning…
and the performance gap widens when doubling the number of bins.
so finally let's examine 2d density estimation for heat maps and contour plots we'll look at mileage versus horsepower in the classic cars data set.
using direct calculation as a baseline we can visualize the difference in estimates when using box and extended box filtering.
we again see that these methods tend to underestimate peaks yet overestimate the size of a distribution.
errors from the deriche method however are imperceptibly low on this color scale.
now if we overlay contours from all three methods we'll find that deriche matches the ground truth…
…while the box filter methods can incur missing contours missed estimated peaks and again other hallucinators.
across bandwidths we see that the deriche method is consistently better.
once again this further improves by using linear binning…
…and increasing the grid size.